Nonparametric graphon estimation

Introduction

Networks are fast becoming part of the modern statistical landscape (Durrett, 2007; Diaconis and Janson, 2008; Bickel and Chen, 2009; Choi, Wolfe and Airoldi, 2012; Fienberg, 2012; Zhao, Levina and Zhu, 2012; Arias-Castro and Grimmett, 2013; Ball, Britton and Sirl, 2013; Choi and Wolfe, 2013). Yet we lack a full understanding of their large-sample properties in all but the simplest settings, hindering the development of models and inference tools that admit theoretical performance guarantees.

In this article we introduce a nonparametric framework for the analysis of networks, which relates to kernel-based random graph models (Janson, 2010; Sussman, Tang and Priebe, 2013), stochastic blockmodels (Airoldi et al., 2008; Rohe, Chatterjee and Yu, 2011), and degree-based models (Chatterjee, Diaconis and Sly, 2011; Bickel, Chen and Levina, 2011). We use this framework to establish consistency of likelihood-based network inference under general conditions, and to show convergence rates across a range of network regimes, from dense to sparse. Our framework thus addresses one of the biggest factors limiting the use of statistical network models in practice: a lack of flexible and transparent analysis tools that admit coherent statistical interpretations (Fienberg, 2012).

Our methodology derives from a large-sample theory tailored to network data, in which well-defined limiting objects play a role akin to the infinite-dimensional functions that underpin classical nonparametric statistics (Bickel and Chen, 2009). An exchangeable stochastic network can be viewed as a partial observation of this limiting object under Bernoulli sampling (Diaconis and Janson, 2008). Hence our theory is closely related to that of generalized linear models (Green and Silverman, 1994) and of contingency tables (Fienberg and Rinaldo, 2012), as well as to nonparametric function approximation. High-dimensional statistical theory in this setting is nascent, and so the linkages we develop below provide for a foundational understanding of nonparametric statistical network analysis.

Model elicitation

A network can be represented by an $n\times n$ data matrix $A$, whose $ij$th entry describes the relation between node $i$ and node $j$ of the network. In the most fundamental setting of graph theory, $A$ is a symmetric, binary-valued contingency table: it is sparse yet structured, with $A_{ij}\in\{0,1\}$ denoting the absence or presence of an edge between nodes $i$ and $j$, and with fixed, structural zeros along the main diagonal.

We call $A$ an adjacency matrix, and model it as a realization of $\binom{n}{2}$ independent Bernoulli trials. Independently for $1\leq i<j\leq n$, we have

Each Bernoulli trial $A_{ij}$ has success probability $p_{ij}$, which in turn we model using a bivariate function termed a graphon that derives from the theory of graph limits (Lovász, 2012).

A graphon is a nonnegative symmetric function, measurable and bounded, that represents a discrete network as an infinite-dimensional analytic object. It is a basic characterization, allowing us to go from the discrete set of probabilities $\{p_{ij}\}_{i<j}$ to a limit object $f\left(x,y\right)$ defined on $(0,1)^{2}$, independently of the network size. Various summaries of the network can be calculated as functionals of the graphon; for example, a network’s degree distribution is characterized by its graphon marginal $\int_{0}^{1}f\left(\cdot,y\right)\,dy$.

To model both dense and sparse networks, we allow the success probabilities $p_{ij}$ appearing in (2.1) to depend on $n$. We link these to a scaled graphon $\rho_{n}f\left(x,y\right)$ through a random sample $\{\xi_{i}\}_{i=1}^{n}$ of uniform variates, via a scale parameter $\rho_{n}>0$ that specifies the expected probability of a network edge:

Observe that $\operatorname{E}A_{ij}=\operatorname{E}_{\xi}p_{ij}=\rho_{n}$ for all $1\leq i<j\leq n$, and so $\rho_{n}$ specifies the sparsity of the generated network. We assume the sequence $\{\rho_{n}\}_{n=2,3,\ldots}$ to be fixed and monotone non-increasing.

This is a canonical model based on exchangeable random networks (Bickel and Chen, 2009; Bickel, Chen and Levina, 2011), and is also strongly related to other statistical modeling paradigms. It relates the infinite-dimensional graphon $f\left(x,y\right)$ to the set of probabilities $\{p_{ij}\}_{i<j}$ sampled via $\xi$. This modeling strategy is similar to time series analysis, where a sampled autocovariance is related to an infinite-dimensional spectral representation. As with an independent increments process, we may think of each $\xi_{i}$ in (2.2) as a latent variable. Furthermore, $\xi_{i}$ is associated with the $i$th network node, acting as a latent random index into the graphon. This reflects the fact that the observed ordering of the network nodes conveys no information.

Similarly, the ordering of a given graphon $f\left(x,y\right)$ along the $x$ and $y$ axes has no inherent meaning; that is, $f\left(x,y\right)$ has a built-in invariance to “rearrangements” of the $x$ and $y$ axes. This is similar to statistical shape analysis, where we seek to describe objects in a manner that is invariant to their orientation in Euclidean space. Thus $f\left(x,y\right)$ represents an equivalence class of all symmetric functions that can be obtained from one another through measure-preserving transformations of $$.

This notion was formalized by Aldous (1981) and Hoover (1979) in the context of exchangeable infinite arrays. Their eponymous theorem asserts that any such array admits a representation in terms of some $f(x,y,\alpha)$. This representation is unique up to measure-preserving transformation (Diaconis and Janson, 2008), and the value of $\alpha$ is not identifiable from a single network observation (Bickel and Chen, 2009). The Aldous–Hoover representation thus relates (2.2) to an exchangeable infinite array $\{A_{ij}\}_{i,j=1}^{\infty}$ of binary random variables, such that for all $n=1,2,\ldots$, all permutations $\Pi$ of $\{1,\ldots,n\}$ and all $a\in\{0,1\}^{n\times n}$, we have that $\Pr(A_{ij}=a_{ij},1\leq i<j\leq n)=\Pr(A_{ij}=a_{\Pi(i)\Pi(j)},1\leq i,j\leq n)$.

By putting an observed $n\times n$ adjacency matrix $A$ in correspondence with a finite set of rows and columns of $\{A_{ij}\}_{i,j=1}^{\infty}$, we arrive at a model for exchangeable networks, or for sub-networks thereof. Exchangeability implies that once we condition on the latent variable $\xi_{i}$ associated to network node $i$, then all linkages $A_{i\cdot}$ to node $i$ are conditionally independent and identically distributed. This follows from de Finetti’s representation of a sum of exchangeable indicator variables (Diaconis, 1977).

Main result

Our main result is that whenever a graphon $f$ is Hölder continuous, and maximum likelihood fitting is used to derive a nonparametric estimator of $f$ from $A$, then this estimator will be consistent as long as \rho_{n}=\omega\bigl{(}n^{-1}\log^{3}n\bigr{)}, and its rate of convergence can be established.

To construct our estimator, we will calculate group averages after forming $k$ groups from $n$ nodes. Any such grouping can be represented as an integer partition of $n$ via a vector $h\in\{2,\ldots,n\}^{k}$, such that $\smash{\sum_{a=1}^{k}}h_{a}=n$. Thus may view $n^{-1}h$ as the probability mass function of a random variable with range $\left\{1,\ldots,k\right\}$, indexed via a cumulative distribution function $H$ and its generalized inverse $H^{-1}$:

The central difficulty in constructing a nonparametric graphon estimator is that we do not know the ordering of our observed adjacency matrix $A$, relative to the ordered sample $\{\xi_{(i)}\}_{i=1}^{n}$ indexing the graphon $f$. We thus define an estimator $\smash{\hat{f}}$ as a composition of two operations: first we re-index the rows and columns of $A$ according to some permutation $\Pi$ of $\{1,\ldots,n\}$, and then we group them in accordance with $H$:

We then define the mean-squared error of $\hat{f}$ relative to $f$ as

where $\mathcal{M}$ is the set of all measure-preserving bijections of the form $\sigma\colon\to$. This error criterion is based on the so-called cut distance in the theory of graph limits (Lovász, 2012), and allows for all possible rearrangements of the axes of $f$ (Choi and Wolfe, 2013).

Any estimator $\hat{f}$ can be viewed as a Riemann sum approximation of $f$, and thus we must understand when such sums converge. Lebesgue’s criterion asserts that a bounded graphon on $(0,1)^{2}$ is Riemann integrable if and only if it is almost everywhere continuous. A sufficient condition is that $f$ is $\alpha$-Hölder continuous for some $0<\alpha\leq 1$, where we write

This assumption ensures that $f$ is uniformly continuous, so that its approximation error can be controlled through Riemann sums.

Under this model specification, we obtain our main result, which we prove in Appendix A.

Assume a sequence of graphon estimators $\hat{f}\left(x,y;h\right)$ is fitted under the model of (2.2), with $k=\omega(1)$ and $\bar{h}=n/k$ the average group size, where

The graphon $f$ is symmetric, bounded away from zero and $\alpha$-Hölder continuous, $0<\alpha\leq 1$;

The scaling sequence $\rho_{n}$ satisfies \rho_{n}=\omega\bigl{(}n^{-1}\log^{3}n\bigr{)}, and $\max_{n}\rho_{n}f$ is bounded away from unity;

Every admissible partition $H$ has group sizes bounded uniformly above and below by $h_{\vee}=o(n)$, $h_{\wedge}=\omega(\log^{1/2}n)$, and may be composed with any permutation $\Pi$ of $\{1,\ldots,n\}$ to yield $\hat{f}\left(x,y;h\right)$.

Suppose furthermore that the minimum effective sample size of every possible fitted grouping, $\binom{h_{\wedge}}{2}\rho_{n}$, and the average effective sample size across all groupings, $\bar{h}^{2}\rho_{n}$, both grow sufficiently rapidly in $n$:

Then if $\hat{f}\left(x,y;h\right)$ is fitted by blockmodel maximum profile likelihood estimation as described in Section 4 below, the mean-squared error of $\hat{f}$ satisfies

The terms appearing in this expression each stem from a different portion of the nonparametric inference problem of graphon estimation, and will be derived and discussed in Section 5–7 below.

Nonparametric graphon approximation via blockmodels

To understand Theorem 3.1, we must first describe how a particular class of statistical network model—the stochastic blockmodel—lends itself naturally to nonparametric approximation. Later, in Section 5, we will establish blockmodel consistency under model misspecification, in settings ranging from dense (Chatterjee, 2012; Choi and Wolfe, 2013) to very sparse networks.

A $k$-community blockmodel $(k,z,\theta)$ is a statistical network model that consists of two main components:

A community assignment function $z\colon\left\{1,\ldots,n\right\}\to\left\{1,\ldots,k\right\}$. This mapping assigns each of $n$ network nodes to exactly one of $k$ groupings or “communities,” each of size $h_{a},1\leq a\leq k$.

A block mean estimator $\theta\colon\left\{1,\ldots,k\right\}^{n}\times\left[0,1\right]^{n\times n}\to\left[0,1\right]^{k\times k}$. This assigns an interaction rate $\theta_{ab}$ to every pair $(a,b)$ of communities, based on the observations $\left\{A_{ij}:i\in z^{-1}(a),j\in z^{-1}(b)\right\}$.

Any community assignment function $z$ thus has two components: a vector $h(z)=\left(h_{1},\ldots,h_{k}\right)$ of community sizes equivalent to some $H$ as defined in (3.1a), and a permutation $\Pi_{z}$ of $\left\{1,\ldots,n\right\}$ that re-orders the set of network nodes prior to applying the quantile function $H^{-1}(\cdot/n)$ as defined in (3.1b). Thus the community to which $z$ assigns node $i$ is determined by the composition $H^{-1}\circ\Pi_{z}$:

Each $z$ thus represents a re-ordering of the network nodes, followed by a partitioning of the unit interval. Each $\theta_{ab}$ in turn describes the expected rate of interaction between the nodes in communities $a$ and $b$.

If $k$ grows with $n$, then the nonparametric properties of blockmodels come to the fore (Rohe, Chatterjee and Yu, 2011; Choi, Wolfe and Airoldi, 2012; Fishkind et al., 2013; Zhao, Levina and Zhu, 2012). In the theory of graph limits (Lovász, 2012), such a model is known as the “blowup” of a weighted graph to the domain $(0,1)^{2}$, or as a “stepfunction approximation” of a given graphon $f\left(x,y\right)$.

There are strong theoretical reasons why an arbitrary graphon should be well approximated by blocks (Lovász, 2012). These reasons stem from a fundamental result in combinatorics known as Szemerédi’s regularity lemma, which cuts across graph theory, analysis and number theory. In our context, this lemma suggests that any sufficiently large graph behaves approximately like a $(k,z,\theta)$-blockmodel for some $k$. However, this value of $k$ may potentially be very large, and so regularizing strategies are needed to infer a blockmodel approximation with good risk properties while requiring relatively few degrees of freedom.

2 Fitting blockmodels to inhomogeneous random graphs

Once $f\left(x,y\right)$ has been specified and a uniform random sample $\{\xi_{i}\}_{i=1}^{n}$ realized, our network reduces to a set of $\binom{n}{2}$ $\operatorname{Bernoulli}(p_{ij})$ trials that are conditionally independent given $\{\xi_{i}\}_{i=1}^{n}$. We refer to this as an inhomogeneous random graph model (Bollobás, Janson and Riordan, 2007) for the observed data matrix $A\in\{0,1\}^{n\times n}$. From (2.2), the conditional log-probability of observing a given adjacency matrix $A$ is

Adopting the notation of Choi, Wolfe and Airoldi (2012), we write the log-likelihood function of a blockmodel $(k,z,\theta)$ with respect to an observed data matrix $A$ as

where $\bar{A}_{ab}$ is the arithmetic average of the values of $A$ in the $(a,b)$th block:

and $h_{a}$ is the size of the $a$th community. Note that this aligns with our earlier definition of $\hat{f}$, and that the quantities $h_{ab}^{2},\bar{A}_{ab},\theta_{ab}$ all depend on the community assignment function $z$. The structural zeros along the main diagonal of $A$ imply that $h_{ab}$ differs for diagonal blocks ($a=b$) relative to off-diagonal blocks. We see from (4.2) that for any fixed assignment $z\in\{1,\ldots,k\}^{n}$, the log-likelihood $L(A;z,\theta)$ of $A$ will be maximized in $\theta\in^{k\times k}$ by taking $\theta_{ab}=\bar{A}_{ab}$. This is because each sample proportion $\bar{A}_{ab}$ is an extended maximum likelihood estimator for its expectation; “extended”, because we include the boundary $\{0,1\}^{k\times k}$ of the parameter space, allowing for the possibility that $\theta_{ab}=\bar{A}_{ab}\in\{0,1\}$. Thus the extended maximum likelihood estimator coincides with the method of moments estimator for $\theta_{ab}$.

Note that (4.2) is a continuous function in $\theta$, and so (by the extreme value theorem) $L(A;z,\theta)$ attains its supremum over the compact set $^{k\times k}$. Thus we “profile out” $\theta$ from the log-likelihood $L(A;z,\theta)$:

Any maximizer of (4.4) over a fixed, non-empty subset $\mathcal{Z}_{k}\subseteq\{1,\ldots,k\}^{n}$ is a maximum profile likelihood estimator (MPLE) of $z$ with respect to $\mathcal{Z}_{k}$. We may equivalently re-cast the problem of likelihood maximization as one of Bernoulli Kullback–Leibler divergence minimization, with

denoting the Kullback–Leibler divergence of a $\operatorname{Bernoulli}(p^{\prime})$ distribution from a $\operatorname{Bernoulli}(p)$ one.

Equipped with this definition, observe that any MPLE $\hat{z}(A,\mathcal{Z}_{k})$ satisfies

Maximizing the profile log-likelihood of (4.4) to obtain an MPLE $\hat{z}(A,\mathcal{Z}_{k})$ is thus equivalent to minimizing the sum of divergences $\sum_{i<j}\operatorname{D}\left(A_{ij}\,\middle|\middle|\,\bar{A}_{z_{i}z_{j}}\right)$. This sum serves as a proxy for its “oracle” counterpart based on the matrix $p\in^{n\times n}$ of Bernoulli parameters of the underlying generative model. This corresponds to an idealized “best blockmodel approximation” of $p$.

With this in mind, we define an “oracle MPLE” $z(p,\mathcal{Z}_{k})$ in direct analogy to (4.5). Let $\bar{p}(z)_{ab}$ denote the arithmetic average of the $h_{ab}^{2}$ elements of $p$ in the $(a,b)$th block induced by $z$:

where we recall that $h_{ab}^{2}$ also depends on the choice of community assignment function $z$. We then have

Observe that neither $\hat{z}(A,\mathcal{Z}_{k})$ nor $\bar{z}(p,\mathcal{Z}_{k})$ is unique, since permuting the community labels $\{1,\ldots,k\}$ does not affect the likelihood of community assignment in (4.5) or (4.7). Even aside from the issue of label switching, we are not guaranteed uniqueness; see Chatterjee, Diaconis and Sly (2011) and Rinaldo, Petrović and Fienberg (2013) for discussion of this issue in the specific context of network modeling, as well as Fienberg and Rinaldo (2012) in the general setting of log-linear models for sparse contingency tables.

Sparse blockmodel consistency under model misspecification

We now establish that an observed matrix $A\in\{0,1\}^{n\times n}$ of binary adjacencies yields “oracle” information on its generative $p\in(0,1)^{n\times n}$ at a rate that depends both on the sparsity of the network and on the speed at which the admissible network community sizes grow with $n$. We show that for suitable sequences of sets $\mathcal{Z}_{k}(n)\subseteq\{1,\ldots,k\}^{n}$ of admissible blockmodels, the maximum profile likelihood assignment method $\hat{z}(A,\mathcal{Z}_{k})$ implies that the likelihood risk of a fitted blockmodel, as measured by summing the divergences $\operatorname{D}\left(p_{ij}\,\middle|\middle|\,\bar{A}_{\hat{z}_{i}\hat{z}_{j}}\right)$, approaches the risk $\sum_{i<j}\operatorname{D}\left(p_{ij}\,\middle|\middle|\,\bar{p}_{z_{i}z_{j}}\right)$ of the best possible blockmodel approximation as $n$ grows large.

Theorem 5.1 (proved in Appendix B) makes this statement precise and provides a set of sufficient conditions, driven primarily by the effective sample size of each fitted block.

For each $n=2,3,\ldots$, let $A\in\{0,1\}^{n\times n}$ be the adjacency matrix of a simple random graph with independent $\operatorname{Bernoulli}(p_{ij})$ edges, and consider a corresponding sequence of $k$-community blockmodel estimators, with $k=k(n)$ a function of $n$. Assume:

The expected edge density $\binom{n}{2}^{-1}\sum_{i<j}p_{ij}(n)$ of $A$ does not approach or $1$ too rapidly in $n$: there exists a monotone non-increasing, strictly positive sequence $\bar{\rho}(n)$, such that for all $n$ sufficiently large, $\bar{\rho}(n)\leq\binom{n}{2}^{-1}\sum_{i<j}p_{ij}(n)\leq 1-\sqrt{\bar{\rho}(n)}$.

Likewise, no block density $\{\bar{p}_{z_{i}z_{j}}(n)\}_{i<j,z\in\mathcal{Z}_{k}(n)}$ approaches or $1$ too rapidly in $n$: there exists a monotone non-increasing, strictly positive sequence $\rho_{\wedge}(n)$, such that $\rho_{\wedge}(n)\leq\bar{\rho}(n)$ and $\rho_{\wedge}(n)\leq\bar{p}_{z_{i}z_{j}}(n)\leq 1-\sqrt{\rho_{\wedge}(n)}$ for all $z\in\mathcal{Z}_{k}(n)$, $1\leq i<j\leq n$ and $n$ sufficiently large.

The sizes $\{h_{z_{i}}(n)\}_{1\leq i\leq n,z\in\mathcal{Z}_{k}(n)}$ of all possible communities grow sufficiently rapidly in $n$: there exists a monotone strictly increasing sequence $h_{\wedge}(n)$ taking values in $\{2,\ldots,\lfloor n/k(n)\rfloor$ such that for all $n$ sufficiently large, $h_{\wedge}(n)\leq\min_{z\in\mathcal{Z}_{k}(n)}\left\{\min_{1\leq i\leq n}h_{z_{i}}(n)\right\}$.

Assume that the sequences $\mathcal{Z}_{k},\bar{\rho},\rho_{\wedge},h_{\wedge}$ are fixed in advance and independent of all other quantities. Let $\bar{h}=n/k\in[1,n]$, and suppose that the minimum effective sample size of every possible fitted block, $\binom{h_{\wedge}}{2}\rho_{\wedge}$, and the average effective sample size across all blocks, $\bar{h}^{2}\bar{\rho}$, both grow sufficiently rapidly in $n$:

Then for all sequences of subsets $\mathcal{Z}_{k}\subseteq\{1,\ldots,k\}^{n}$ that respect condition 3, we have as $n\rightarrow\infty$ that for any choice of $z\in\mathcal{Z}_{k}$, deterministic or random,

For $\hat{z}(A,\mathcal{Z}_{k})=\operatorname{argmax}_{z\in\mathcal{Z}_{k}}\sum_{i<j}\left\{A_{ij}\log\bar{A}_{z_{i}z_{j}}+\left(1-A_{ij}\right)\log\left(1-\bar{A}_{z_{i}z_{j}}\right)\right\}$,

These results also hold marginally with respect to the model of (2.2).

Theorem 5.1 is significant because it gives conditions under which the excess risk of a fitted blockmodel converges to zero, implying that blockmodel parameters can be estimated consistently even when the true generative model giving rise to $A$ is unknown. It predicts different rates of convergence for different network sparsity regimes. Depending on the growth of $k$ with $n$, either the first or the second of two rate terms in (5.2) will dominate.

We may summarize these regimes as follows:

Dense networks: If $\rho_{\wedge}$ and $\bar{\rho}$ remain constant in $n$, and $k$ grows with $n$ as $k=\mathcal{O}(n^{3/4})$, then Theorem 5.1 predicts a convergence rate of at least $\sqrt{\log(n)/n}$. If instead $k$ grows like $n^{\delta}$ for $3/4<\delta<1$, then this rate will decrease to $\log n/n^{2(1-\delta)}$.

Sparse networks: If $\rho_{\wedge}$ and $\bar{\rho}$ decrease like $n^{-2\gamma}$ for $0<\gamma<1/2$, and $k=\mathcal{O}(n^{3/4-\gamma/2})$, then Theorem 5.1 predicts the rate $\log(n)^{3/2}/n^{1/2-\gamma}$. If $k$ grows like $n^{\delta}$ for $3/4-\gamma/2<\delta<1-\gamma$, then this rate will decrease to $\log n/n^{2(1-\delta-\gamma)}$.

Ultra-sparse networks: If $\rho_{\wedge}$ and $\bar{\rho}$ decrease like $\log(n)^{3+\beta}/n$ for $\beta>0$, then Theorem 5.1 predicts rate $\log(n)^{-\beta/2}$ whenever $k=\mathcal{O}(n^{1/2})$, matching the regime of Choi, Wolfe and Airoldi (2012).

In each of these cases, the given conditions on $\rho_{\wedge}$ can be relaxed accordingly.

Theorem 5.1 is the first such result known for sparse or ultra-sparse networks—those for which $\bar{\rho}=o(1)$, so that the average number of connections per node can grow sublinearly, here as slowly as logarithmically in $n$. This complements the recent result of Choi and Wolfe (2013) for fixed-$k$ fitting of dense bipartite graphs—those for which $\rho_{\wedge}$ and $\bar{\rho}$ remain constant, so that the average number of connections per node grows linearly in $n$. Theorem 5.1 extends this regime, allowing for the growth of $k$ with $n$, while also yielding an improved convergence rate of $\sqrt{\log(k)/n}$ for dense graphs.

To understand why Theorem 5.1 holds in this setting, we begin by conditioning on a choice of community assignment function $z$. Blocks of network edges then comprise independent sets of independent Bernoulli trials. Conditionally upon $z$, sample proportions $\bar{A}_{z_{i}z_{j}}\,|\,z$ of these blocks are thus independent Poisson–Binomial variates. Without additional restrictions, however, a fitted block could be any size—even as small as a single Bernoulli trial. Thus it is necessary to constrain the set $\mathcal{Z}_{k}\subseteq\{1,\ldots,k\}^{n}$ of admissible blockmodels, and also to constrain the allowable global and local sparsity of the network, so that the effective sample size of every possible $\bar{A}_{z_{i}z_{j}}\,|\,z$ grows in $n$. This ensures that all block-wise sample proportions $\bar{A}_{z_{i}z_{j}}\,|\,z$ behave like Normal variates in the large-sample limit, when appropriately standardized.

There are then two main technical challenges:

Double randomness: While every $\bar{A}_{z_{i}z_{j}}\,|\,z$ is amenable to analysis, choosing $\hat{z}$ by profile likelihood maximization introduces “double randomness,” coupling all blocks and precluding a direct analysis of $\bar{A}_{\hat{z}_{i}\hat{z}_{j}}$. Instead, we take the approach of Choi, Wolfe and Airoldi (2012), and show that results for $\bar{A}_{z_{i}z_{j}}\,|\,z$ hold uniformly for any choice of $z$ — and therefore that they also hold for $\bar{A}_{\hat{z}_{i}\hat{z}_{j}}$.

Likelihood zeros: The assumption that all $p_{ij}\in(0,1)$ ensures that each $\operatorname{D}\left(p_{ij}\,\middle|\middle|\,\bar{p}_{z_{i}z_{j}}\right)$ is finite. However, $\operatorname{D}\left(p_{ij}\,\middle|\middle|\,\bar{A}_{\hat{z}_{i}\hat{z}_{j}}\right)$ will fail to be finite if $\bar{A}_{\hat{z}_{i}\hat{z}_{j}}\in\{0,1\}$, in which case the $(\hat{z}_{i},\hat{z}_{j})$th block has saturated. Such blocks add to the likelihood; their parameters are not estimable (Fienberg and Rinaldo, 2012). The theorem conditions allow us to control the probability of these likelihood zeros, by requiring the effective sample size of each block to grow sufficiently rapidly in $n$.

This latter point is particularly important, since only values in the interior of the parameter space $^{k\times k}$ are estimable (Fienberg and Rinaldo, 2012, Theorem 7). As in the case of additional structural zeros (Fienberg and Rinaldo, 2012, Corollary 8), the Fisher information matrix will be rank-deficient, and the degrees of freedom must be adjusted accordingly in order to obtain correct inferential conclusions. This explains why the random denominator term is necessary in the left-hand side of (5.2).

We may connect this understanding to the three sparsity regimes described above: the case of dense networks, corresponding to the setting of exchangeable random graphs; that of sparse networks, where the density of network edges $\smash{\binom{n}{2}^{-1}}\sum_{i<j}p_{ij}$ decays as some power of $n$; and that of ultra-sparse networks, where the edge density decays at a rate approaching $\log(n)/n$. This is the so-called connectivity threshold, above which an inhomogeneous random graph will be fully connected with probability approaching $1$ as $n\rightarrow\infty$ (Alon, 1995). If the edge density were instead to decay at a rate of $1/n$—the extremely sparse setting of Bollobás and Riordan (2009)—then the resulting networks would fail in general to be connected, and Poisson rather than Normal limiting behavior would hold for each block (Olhede and Wolfe, 2013).

From blockmodels to smooth graphon estimation

We now present our final result leading to consistent graphon estimation. To go beyond conditional estimation of inhomogeneous random graphs via blockmodels, we will assume additional structure via graphon smoothness. This smoothness will in turn allow us to control estimation risk, by sending the main term in Theorem 5.1 to zero.

A blockmodel first orders the rows and columns of $A$, and then groups its entries according to a vector of community sizes $h\in\{2,\ldots,n\}^{k}$. This specifies a partition $H$ in accordance with (3.1a), which in turn induces a piecewise-constant approximation of the graphon $f\left(x,y\right)$ along blocks. To see this, define the domain $\omega_{ab}\subseteq[0,1)^{2}$ of the $(a,b)$th block as

If $f\left(x,y\right)$ is smooth as well as bounded, then results from approximation theory allow the error $\|f-\bar{f}\|$ to be controlled in any $L_{p}$ norm, as a function of the maximum over all block diameters $(h_{a}^{2}+h_{b}^{2})^{1/2}/n$ for $1\leq a,b\leq k$ (DeVore, 1998, see also Lemma C.6).

Recall from (4.1) that any blockmodel community assignment vector $z$ is a composition $H^{-1}\circ\Pi_{z}$ for some partition $H$ of $$and permutation$\Pi_{z}$of$\{1,\ldots,n\}$, so that$z_{i}=H^{-1}\left\{\Pi_{z}(i)/n\right\},1\leq i\leq n$. From (4.6), we may express$\bar{p}(z)$for any$1\leq a,b\leq k$ as

Thus $\bar{p}(z)_{ab}$ is an average over $h_{ab}^{2}$ graphon evaluations f\,\bigl{(}\xi_{\Pi_{z}^{-1}(i)},\xi_{\Pi_{z}^{-1}(j)}\bigr{)}, since the model of (2.2) asserts that $p_{ij}(n)\propto f\left(\xi_{i},\xi_{j}\right)$. These evaluations occur at random points determined by $\{\xi_{1},\ldots\xi_{n}\}$ according to the inverse of the permutation $\Pi_{z}$, while $H$ determines the size of each block.

From this simple observation, we will show that it is possible to relate $\bar{p}(z)_{ab}$ to $f\left(x,y\right)$ by choosing an “oracle” permutation $\Pi_{z}(i)$ whose inverse yields the ordered sample $\{\xi_{(1)},\ldots\xi_{(n)}\}$. To see this, first note that whenever the Hölder condition of (3.2) is satisfied, we have by Lemma C.7 that

because each $\xi_{(i)}$ converges in probability to its expectation $i/\left(n+1\right)$ at a rate no worse than $n^{-1/2}$, and (3.2) relates this to \bigl{|}f\left(\xi_{(i)},\xi_{(j)}\right)-f\,\bigl{(}\frac{i}{n+1},\frac{j}{n+1}\bigr{)}\bigr{|}. Now take $\Pi_{z}(i)=(i)^{-1}$, where $(i)^{-1}$ denotes the rank of $\xi_{i}$ from smallest to largest, and observe that f\,\bigl{(}\xi_{\Pi_{z}^{-1}(i)},\xi_{\Pi_{z}^{-1}(j)}\bigr{)} evaluates to $f\left(\xi_{(i)},\xi_{(j)}\right)$.

The key point is that when $f$ is $\alpha$-Hölder continuous, then convergence of the ordered sample $\{\xi_{(i)}\}_{i=1}^{n}$ governs convergence of the random averages comprising $\bar{p}(z)_{ab}$ in (6.2). Indeed, if $h_{\vee}$ uniformly upper-bounds the largest possible community size, then by Lemma C.5, we have that

where we recall from (6.1) that ${\bar{f}}\left(x,y;h\right)$ is the local block average of $f$.

As a consequence, we can control the oracle estimation risk featured in Theorem 5.1 as follows.

Assume in the scaled exchangeable graph model of (2.2) that:

The graphon $f$ is a positive, symmetric function on $(0,1)^{2}$, and is $\alpha$-Hölder continuous, $0<\alpha\leq 1$;

Furthermore, $f$ is bounded away from zero and $\max_{n}\rho_{n}f$ is bounded away from unity;

Each set $\mathcal{Z}_{k}(n)\subseteq\{1,\ldots,k\}^{n}$ of admissible blockmodel assignments has the following property: If $H$ is generated by some $z\in\mathcal{Z}_{k}$, then $H^{-1}\circ\Pi\in\mathcal{Z}_{k}$ for every permutation $\Pi$ of $\{1,\ldots,n\}$.

Then for $h_{\vee}(n)$ the largest community size in each $\mathcal{Z}_{k}(n)$, the oracle likelihood risk in Theorem 5.1 satisfies

We prove this theorem in Appendix C by using the oracle choice of permutation $(\cdot)^{-1}$ to upper-bound the risk via a block approximation ${\bar{f}}\left(x,y;h\right)$ of $f\left(x,y\right)$, based on some $z^{*}$ which achieves the minimum in (6.3). Conditions 1 and 2 are then sufficient to guarantee the claimed rate of approximation. Condition 3 ensures that $H^{-1}\circ(\cdot)^{-1}\in\mathcal{Z}_{k}$, since we do not know $z^{*}$ or the requisite ordering $(\cdot)^{-1}$ in advance.

Rates of convergence

We see directly that the rate of convergence in Theorem 6.1 depends on the Hölder continuity of $f$ in two ways: through the convergence of the ordered sample $\{\xi_{(i)}\}_{i=1}^{n}$ (variance), and through the rate at which $h_{\vee}/n$ goes to zero in $n$ (bias). This rate is also self-scaling relative to the sparsity of the network, as it does not depend on $\rho_{n}$.

In contrast, Theorem 5.1 depends strongly both on the network sparsity factor $\rho_{n}$, as well as the minimum and average admissible block sizes, $h_{\wedge}$ and $\bar{h}$. The conditions of Theorem 5.1 ensure that excess blockmodel risk can be controlled under model misspecification, enabling groupings of nodes with good risk properties to be estimated, despite the variability of the data.

Together, the results of Theorems 5.1 and 6.1 enable us to establish mean-square graphon consistency at the rates indicated in Theorem 3.1, namely

The first two terms come directly from Theorem 5.1, while the third is from Theorem 6.1. The final term comes from relating the discrete quantities featured in these theorems to the graphon itself, and is driven in part by the fact that we do not know the ordering of the data relative to the $\operatorname{Uniform}(0,1)$ variates $\{\xi_{i}\}_{i=1}^{n}$ by which the graphon is sampled. The \mathcal{O}\bigl{(}n^{-1/2}\bigr{)} variance of the ordered sample $\{\xi_{(i)}\}_{i=1}^{n}$ subsequently appears, and is modulated by the regularity of the graphon through its Hölder continuity exponent $\alpha$.

Conclusion

In this article we have established a number of new results within a nonparametric framework for network inference, based on graphons as natural limiting objects. Understanding graphons as analytic objects, as well as the behavior of dense and sparse networks based on them, is fundamental to advancing our nonparametric understanding of networks.

To this end, we have established consistency of graphon estimation under general conditions, giving rates which include the important practical setting of sparse networks. By treating dense and sparse stochastic blockmodels with a growing number of classes, under model misspecification, our results improve substantially upon what is currently known in the literature.

Our results link strongly to approximation theory, nonparametric function estimation, and the theory of graph limits, and thus provide for a foundational understanding of nonparametric statistical network analysis.

Appendix A Proof of Theorem 3.1 and its lemmas

We note from Lemma A.1 that for $(x,y)\in(0,1)^{2}$

Recalling the definition of $\bar{A}_{ab}$, we see that uniformly for all choices of $H$ and $\Pi$, and for all $1\leq a,b\leq k$, we have $0\leq\operatorname{E}\bar{A}_{ab}\leq\rho_{n}\sup_{(x,y)\in(0,1)^{2}}f\left(x,y\right)$ and $0\leq\operatorname{E}\bar{A}_{ab}^{2}\leq\rho_{n}^{2}\sup_{(x,y)\in(0,1)^{2}}f^{2}\left(x,y\right)$.

Since $f$ is by hypothesis Hölder continuous on a bounded domain, it is bounded, and thus \bar{A}_{ab}={\cal O}_{P}\bigl{(}\rho_{n}\bigr{)} and \bar{A}_{ab}^{2}={\cal O}_{P}\bigl{(}\rho_{n}^{2}\bigr{)} by Markov’s inequality. We will thus expand the squared error term in the integrand of the graphon mean-squared error pointwise, using the fact that the error term should be evaluated at the infimum over measure preserving bijections. Therefore this error be upper-bounded by its evaluation at some $\sigma^{\ast}\in\mathcal{M}$, which we will choose in accordance with the proof of Lemma A.3 below:

where the last two lines follow from Lemmas A.2 and C.9, respectively. By Lemma A.3, we have

uniformly in $z$. The conditions of Theorem 3.1 are sufficient for Theorems 5.1 and 6.1 to hold, and so if $\hat{f}$ is fitted by maximum profile likelihood, then we may substitute terms from Theorems 5.1 and 6.1 to obtain

A.2 Auxiliary lemmas needed for Theorem 3.1

Assume the setting of Theorem 3.1. Then $\operatorname{E}\hat{\rho}_{n}=\rho_{n}$, \operatorname{var}\hat{\rho}_{n}={\cal O}\bigl{(}\rho_{n}^{2}/n\bigr{)}.

The necessary marginal variances and covariances can then be established hierarchically:

Since $\operatorname{E}f\left(\xi_{i},\xi_{j}\right)f\left(\xi_{k},\xi_{l}\right)=\iint_{(0,1)^{2}}f^{2}(x,y)\,dx\,dy$ if $i=k$ and $j=l$, and \bigl{\{}\iint_{(0,1)^{2}}f(x,y)\,dx\,dy\bigr{\}}^{2} if $i\neq k$ and $j\neq l$, we obtain when either $i\neq k$ or $j\neq l$ that

The order term of ${\cal O}(\rho_{n}^{2}/n)$ follows, as $\rho_{n}^{2}/n\geq\rho_{n}/n^{2}\Leftrightarrow\rho_{n}\geq 1/n$, since \rho_{n}=\omega\bigl{(}n^{-1}\log^{3}n\bigr{)}. ∎

Assume the setting of Theorem 3.1. Then for any $z\in{\cal Z}_{k}$,

We first treat the numerator of (A.1), whose infimum is over $\mathcal{M}$, the set of all measure-preserving bijective maps of the form $\sigma\colon\to$. We may write

since $\hat{f}$ is constant on blocks. Observe that for each individual summand in (A.2), we may write

We now restrict our choice of $\sigma\in\mathcal{M}$ to satisfy the following property:

for some permutation $\Pi$ of $\{1,\ldots,n\}$. Such a choice of measure-preserving bijection can always be made, as it simply partitions the unit interval into $n+1$ subintervals of the form $\left[(i-1)/n,i/n\right),1\leq i\leq n$, and permutes their order in accordance with $\Pi$. We make this choice in order to preserve the Hölder continuity of $f$ on each domain $(x,y)\in\left(\frac{i-1}{n},\frac{i}{n}\right)\times\left(\frac{j-1}{n},\frac{j}{n}\right)$, as will be shown below.

Thus we may write, combining (A.2)–(A.6),

with $S_{n}$ the set of permutations of $\{1,\ldots,n\}$. From Lemma A.4 we then obtain

where $\xi_{\left(\Pi\left\{i\right\}\right)}$ is the $\Pi(i)$th element of the ordered sample $\{\xi_{(i)}\}_{i=1}^{n}$. Starting from (A.5), we then have

where we have chosen $\Pi=\left(\cdot\right)^{-1}\circ\Pi_{z}^{-1}$, so that $\Pi(i)=\left(\Pi_{z}^{-1}\{i\}\right)^{-1}$, with $\left(i\right)^{-1}$ the rank of $\xi_{i}$, from smallest to largest. This choice allows us to match each $\xi_{\left(\Pi\left\{i\right\}\right)}$ to the corresponding group assignment $z_{i}$ of the $i$th network node. To see this, recall from (4.1) that $z_{i}=H^{-1}\left\{\Pi_{z}(i)/n\right\},1\leq i\leq n$, and from (4.3) and (C.6) respectively that

Note that $\bar{p}(z)_{ab}=\operatorname{E}\left\{\bar{A}(z)_{ab}\,|\,\xi,z\right\}$. Thus we relate each $p_{ij}=\rho_{n}f\left(\xi_{i},\xi_{j}\right)$ to the average $\bar{A}(z)_{z_{i}z_{j}}$ of the block to which it is assigned by $z$.

Continuing from (A.6), we appeal to Lemma A.5 to bound the diagonal term, thereby obtaining

Lemma A.6 yields the denominator of (A.1), and the result follows by taking the ratio of these terms. ∎

Assume the setting of Theorem 3.1. Then for $1\leq i,j\leq n,\quad(a,b):\bar{A}(z)_{ab}\notin\{0,1\}$

The result follows from a Taylor series of the integrand of (A.7), which we will show to converge everywhere on the domain of integration, as long as $\bar{A}(z)_{ab}\notin\{0,1\}$. We begin by noting that whenever $f\in\operatorname{\textrm{H\"{o}lder}}^{\alpha}(M)$, we have from Lemma C.7 that for all $(x,y)\in\left(\frac{i-1}{n},\frac{i}{n}\right)\times\left(\frac{j-1}{n},\frac{j}{n}\right)$,

From Markov’s inequality, f\left(\xi_{(i)},\xi_{(j)}\right)=f\left(x,y\right)+{\cal O}_{P}\bigl{(}n^{-\alpha/2}\bigr{)} for every fixed $(x,y)$ in the domain of interest. Thus the following Taylor series holds whenever $f\in\operatorname{\textrm{H\"{o}lder}}^{\alpha}(M)$ and $\bar{A}(z)_{ab}\notin\{0,1\}$:

To bound the second term in (A.8), let $l=\inf_{x\in(0,1)}f(x,x)$ and $u=\sup_{x\in(0,1)}f(x,x)$. Since $\bar{A}(z)_{aa}\notin\{0,1\}$, we may bound the magnitudes of $\log\bar{A}(z)_{aa},\log\left\{1-\bar{A}(z)_{aa}\right\}$ via $\log\binom{h_{a}}{2}\leq\log\binom{h_{\vee}}{2}$. Then

The first two terms in (A.9) are bounded by hypothesis, and then we apply Markov’s inequality to (A.8). ∎

Let $l=\inf_{x\in(0,1)}f(x,x)$ and $u=\sup_{x\in(0,1)}f(x,x)$. Since $\bar{A}(z)_{aa}\notin\{0,1\}$, we may bound the magnitudes of $\log\bar{A}(z)_{aa}$ and $\log\left\{1-\bar{A}(z)_{aa}\right\}$ via $\log\binom{h_{a}}{2}\leq\log\binom{h_{\vee}}{2}$. We bound the expectation of each summand in (A.10) for $1\leq i\leq n$

The result then follows from linearity of expectation and Markov’s inequality, as per Lemma A.4. ∎

We start by discretizing the integral. We therefore write that

where the latter term may be bounded using the technique of Lemma A.4, yielding

Note $\sum_{i,j:\bar{A}(z)_{z_{i}z_{j}}\notin\{0,1\}}p_{ij}=2\sum_{i<j:\bar{A}(z)_{z_{i}z_{j}}\notin\{0,1\}}p_{ij}+\sum_{1\leq i\leq n:\bar{A}(z)_{z_{i}z_{j}}\notin\{0,1\}}p_{ii}$, so that

Applying Markov’s theorem and combining the result with (A.11) then yields the stated result. ∎

Appendix B Proof of Theorem 5.1 and lemmas

The proof is divided into four steps, with each the subject of a technical lemma proved in Section B.2.

Lemma B.1 yields the key first step, which is to relate $\operatorname{D}\left(p_{ij}\,\middle|\middle|\,\bar{A}_{z_{i}z_{j}}\right)$ to $\operatorname{D}\left(p_{ij}\,\middle|\middle|\,\bar{p}_{z_{i}z_{j}}\right)$ for any $z\in\mathcal{Z}_{k}$, assuming that $\bar{A}_{z_{i}z_{j}}\notin\{0,1\}$. This ensures that both terms are finite, and hence comparable. To obtain sufficient variance reduction in this setting, every $\bar{A}_{z_{i}z_{j}}$ must concentrate to its mean $\bar{p}_{z_{i}z_{j}}$, in that the ratio of mean to standard deviation must shrink. The minimum effective block sample size $\binom{h_{\wedge}}{2}\rho_{\wedge}$ must grow quickly enough that this takes place, even for the sparsest of all possible fitted blocks.

Assume conditions 1–3 of Theorem 5.1, and that \binom{h_{\wedge}}{2}\rho_{\wedge}=\omega\bigl{(}\log\binom{h_{\wedge}}{2}\bigr{)}. Then

Our next step relies on controlling $\Pr(\bar{A}_{z_{i}z_{j}}\!\in\!\{0,1\})$ uniformly in $z$, via Lemma B.2.

Assume conditions 1–3 of Theorem 5.1. Then

This result shows that the set of all $\bar{A}_{z_{i}z_{j}}\in\{0,1\}$ has vanishing relative cardinality relative to $\sum_{i<j}p_{ij}$, no matter which $z\in\mathcal{Z}_{k}$ is chosen. It is a direct consequence of condition 3 of Theorem 5.1, which ensures that the minimum fitted block size is uniformly lower-bounded by $h_{\wedge}=\omega(1)$.

Lemma B.2 has two immediate consequences. First, we may apply it to conclude that

Second, it enables us to substitute for the term $\sum_{i<j:\bar{A}_{z_{i}z_{j}}\notin\{0,1\}}\operatorname{D}\left(p_{ij}\,\middle|\middle|\,\bar{p}_{z_{i}z_{j}}\right)$ in Lemma B.1 as follows.

Assume conditions 1–3 of Theorem 5.1. Then uniformly for all $z\in\mathcal{Z}_{k}$,

Thus whenever all of the above quantities are $o_{P}(1)$, we may combine Lemmas B.1 and B.3 with (B.1) to obtain our first claimed result: for any choice of $z\in\mathcal{Z}_{k}$, deterministic or random, we have that

whenever conditions 1–3 of Theorem 5.1 hold, \binom{h_{\wedge}}{2}\rho_{\wedge}=\omega\bigl{(}\log\binom{h_{\wedge}}{2}\bigr{)} and the argument of the right-hand side of (B.2) is $o_{P}(1)$. Under these conditions, the numerator term of (B.2), when scaled by $\sum_{i<j}p_{ij}$, converges in probability to and hence in law, whereas (B.1) converges in probability to a non-zero constant. Thus by Slutsky’s theorem, their ratio converges in law, and hence also in probability as per (B.2). Separating terms on the left-hand side of (B.2), and then multiplying the latter numerator term by $\sum_{i<j}p_{ij}/\sum_{i<j}p_{ij}$, we obtain the first result of result of Theorem 5.1, as stated in (5.1).

We now establish sufficient conditions for (B.2). We see immediately that \binom{h_{\wedge}}{2}\rho_{\wedge}=\omega\bigl{(}\log(1/\bar{\rho})\bigr{)} must hold. Since Lemma B.1 requires that \binom{h_{\wedge}}{2}\rho_{\wedge}=\omega\bigl{(}\log\binom{h_{\wedge}}{2}\bigr{)}, we obtain the combined requirement

To see that this condition will be satisfied if the effective sample size of every possible fitted block is \omega\bigl{(}\log n\bigr{)}, first note that $h_{\wedge}\leq n$, and so \log h_{\wedge}^{2}=\mathcal{O}\bigl{(}\log n\bigr{)}. Now observe that because $\rho_{\wedge}\leq\bar{\rho}$, it follows that h_{\wedge}^{2}\rho_{\wedge}=\omega\bigl{(}\log h_{\wedge}^{2}\bigr{)} implies h_{\wedge}^{2}\bar{\rho}=\omega\bigl{(}\log h_{\wedge}^{2}\bigr{)}, or equivalently, \log(1/\bar{\rho})=o\bigl{(}\log(h_{\wedge}^{2}/\log h_{\wedge}^{2})\bigr{)}. Since $h_{\wedge}\leq n$, this in turn implies \log(1/\bar{\rho})=o\bigl{(}\log n\bigr{)}. Thus h_{\wedge}^{2}\rho_{\wedge}=\omega\bigl{(}\log n\bigr{)} implies (B.3) as claimed.

To achieve convergence in probability, (B.2) also requires n^{2}\bar{\rho}=\omega\bigl{(}\log\left|\mathcal{Z}_{k}\right|+\binom{k+1}{2}\bigr{)}. To simplify this requirement and obtain a sufficient condition, observe that $\log\left|\mathcal{Z}_{k}\right|\leq n\log k$, since $\mathcal{Z}_{k}\subseteq\{1,\ldots,k\}^{n}$. Now write $\binom{k+1}{2}=k^{2}\left\{1/2+\mathcal{O}(1)\right\}$, and let $\bar{h}=n/k$. From these simplifications we obtain \bar{\rho}=\omega\bigl{(}\log(n/\bar{h})/n+\bar{h}^{-2}\bigr{)}, which is implied by \bar{h}^{2}\bar{\rho}=\omega\bigl{(}\max\left\{\bar{h}^{2}/n,1\right\}\log n\bigr{)}.

Finally, observe that since the results above hold uniformly over all $z\in\mathcal{Z}_{k}$, they also hold for $z=\hat{z}(A,\mathcal{Z}_{k})$, the maximum profile likelihood estimator of $z$. The following lemma relates this choice to its oracle counterpart $\bar{z}(p,\mathcal{Z}_{k})$—the best choice of $z\in\mathcal{Z}_{k}$—enabling us to strengthen (B.2).

Assume conditions 1 and 2 of Theorem 5.1. Then it follows from the arguments of Theorems 2 and 3 of Choi, Wolfe and Airoldi (2012) that for any $\hat{z}(A,\mathcal{Z}_{k})$ and $\bar{z}(p,\mathcal{Z}_{k})$ as per (4.5) and (4.7),

Since $\bar{z}(p,\mathcal{Z}_{k})$ results in the minimum value of $\sum_{i<j}\operatorname{D}\left(p_{ij}\,\middle|\middle|\,\bar{p}_{z_{i}z_{j}}\right)$, this difference is nonnegative. Its convergence in probability to when suitably normalized is due to the maximizing properties of $\hat{z}(A,\mathcal{Z}_{k})$ and $\bar{z}(p,\mathcal{Z}_{k})$. Thus we conclude that $\hat{z}(A,\mathcal{Z}_{k})$ serves as an empirical proxy for $\bar{z}(p,\mathcal{Z}_{k})$.

To complete the proof, set $z=\hat{z}(A,\mathcal{Z}_{k})$ in (B.2) and combine it with Lemma B.4. Comparing terms, we see that the latter’s will dominate the rate of convergence, and so we upper-bound them using $\bar{h}=n/k=\omega(1)$, subadditivity of the square root and the fact that $\binom{n}{2}/\binom{k+1}{2}\leq\bar{h}^{2}$. We thus obtain

where the final line follows because $\log(n/\bar{h})=o(n\bar{\rho})$ is needed for (B.4) to be $o_{P}(1)$, whereas $\rho_{\wedge}\leq\rho<1/2$ implies that \log(1/\rho_{\wedge})^{2}>\log(2)^{2}=\omega\bigl{(}\,\log(n/\bar{h})/(n\bar{\rho})\,\bigr{)}. Thus we have derived the claimed rate of convergence, with a sufficient condition being that \bar{h}^{2}\bar{\rho}=\omega\bigl{(}\max\left\{\bar{h}^{2}/n,1\right\}\log^{3}n\bigr{)}, since together \bar{h}^{2}\bar{\rho}=\omega\bigl{(}\log n\bigr{)} and \rho=\omega\bigl{(}\log(n)^{3}/n\bigr{)} imply that (B.4) is $o_{P}(1)$.

To complete the proof of Theorem 5.1, we now re-interpret the above results under the scaled exchangeable random graph model of (2.2). Lemmas B.1–B.4 then hold for every realized value of $\xi$, and thus the implicit conditioning on $\xi$ inherent to these results can be removed. Specifically, in Lemmas B.1 and B.4, we may marginalize (B.7) and (B.12) respectively via the law of total probability, noting that their right-hand sides do not depend on $\xi$. For Lemmas B.2 and B.3, we simply note that the bound of (B.8) holds for all $\xi$.

B.2 Proofs and auxiliary lemmas needed for Theorem 5.1

Since (B.5) is a sum of Kullback–Leibler divergences, it is nonnegative. To show its convergence when suitably normalized, we appeal to Lemma B.5 below, which implies the following under conditions 1–3 of Theorem 5.1 and the hypothesis \binom{h_{\wedge}}{2}\rho_{\wedge}=\omega\bigl{(}\log\binom{h_{\wedge}}{2}\bigr{)}:

For every $\epsilon>0$, eventually in $n$ and with $1^{+}/2$ approaching arbitrarily closely to $1/2$,

where (B.6) follows as $\epsilon\sum_{i<j}p_{ij}\geq 0$ and $(1^{+}/2)\binom{k+1}{2}\geq 0$ eventually in $n$, and (B.7) follows from condition 1 of Theorem 5.1, by which $\sum_{i<j}p_{ij}(n)\geq\binom{n}{2}\,\bar{\rho}(n)$ eventually in $n$. ∎

We will bound $\Pr(\bar{A}_{z_{i}z_{j}}\in\{0,1\})$ uniformly in $z$. Observe that for any $1\leq a\leq b\leq k$, conditionally on any $z\in\mathcal{Z}_{k}$, we have by the arithmetic–geometric mean inequality that

Conditions 2 and 3 of Theorem 5.1 stipulate that for every pair $(a,b)$ and every $z\in\mathcal{Z}_{k}$, eventually in $n$, $\rho_{\wedge}(n)\leq\bar{p}_{ab}(n)\leq 1-\sqrt{\rho_{\wedge}(n)}$ and $h_{\wedge}(n)\leq h_{a}(n)$. Hence (B.8) implies that, eventually in $n$, for $1\leq a\leq b\leq k$

Since the conditional probability \Pr\big{(}\bar{A}_{z_{i}z_{j}}\in\{0,1\}\,|\,Z=z\big{)} is upper-bounded by (B.9) uniformly for every value of $z\in\mathcal{Z}_{k}$, this same bound also holds after marginalizing out $Z$. Thus, eventually in $n$,

Applying Markov’s inequality, we see that for any $\epsilon>0$, eventually in $n$,

where the second inequality follows directly from (B.10), the third inequality follows from condition 1 of Theorem 5.1, by which $\sum_{i<j}p_{ij}(n)\geq\binom{n}{2}\,\bar{\rho}(n)$ eventually in $n$, and the final inequality follows from the fact that \log\bigl{\{}(1-\rho_{\wedge})^{\binom{h_{\wedge}}{2}}\bigr{\}}=\binom{h_{\wedge}}{2}\log(1-\rho_{\wedge})\leq-\binom{h_{\wedge}}{2}\rho_{\wedge}. ∎

First, we express the term of interest as a sum of nonnegative random variables:

To show the claimed convergence in probability, we write

In the notation of Choi, Wolfe and Airoldi (2012), define for any fixed $z\in\mathcal{Z}_{k}$

where the implication follows directly from the definition of the “oracle” MPLE in $\bar{z}(p,\mathcal{Z}_{k})$ in (4.7). Thus

By construction, since $\bar{z}(p,\mathcal{Z}_{k})$ maximizes $\bar{L}(z)$ over $\mathcal{Z}_{k}$, this difference is nonnegative. Similarly, from (4.5) we see that $\hat{z}(A,\mathcal{Z}_{k})$ maximizes $L(A;z)$ over $\mathcal{Z}_{k}$, and so $L(A;\hat{z})-L(A;\bar{z})\geq 0$. Hence,

and so the result will follow from (B.11) if we can show that $\left|\bar{L}(\bar{z})-L(A;\bar{z})\right|$ and $\left|L(A;\hat{z})-\bar{L}(\hat{z})\right|$ both converge in probability to zero when suitably renormalized. We accomplish this in the manner of Choi, Wolfe and Airoldi (2012, Theorem 2), who establish that $\max_{z\in\mathcal{Z}_{k}}\left|\bar{L}(z)-L(A;z)\right|/\sum_{i<j}p_{ij}$ converges as required. Since this result holds for the maximum over all $z\in\mathcal{Z}_{k}$, then it must also hold for both $\hat{z}$ and $\bar{z}$, and we can therefore apply this same result twice.