Co-clustering separately exchangeable network data

Introduction

Blockmodels are popular tools for network modeling that see wide and rapidly growing use in analyzing social, economic and biological systems; see Zhao, Levina and Zhu (2011) and Fienberg (2012) for recent overviews. A blockmodel dictates that the probability of connection between any two network nodes is determined only by their respective block memberships, parameterized by a latent categorical variable at each node.

Fitting a blockmodel to a binary network adjacency matrix yields a clustering of network nodes, based on their shared proclivities for forming connections. More generally, fitting a blockmodel to any binary array involves partitioning it into blocks. In this way, blockmodels represent a piecewise-constant approximation to a latent function that generates network connection probabilities. This in turn can be viewed as a histogram-like approximation to a nonparametric generative process for binary arrays; fitting such models is termed co-clustering [Flynn and Perry (2012), Rohe and Yu (2012)].

This article analyzes the performance of stochastic blockmodels for co-clustering under model misspecification, assuming only an underlying generative process that satisfies the condition of separate exchangeability [Diaconis and Janson (2008)]. This significantly generalizes known results for the blockmodel and its co-clustering variant, which have been established only recently under the requirement of correct model specification [Bickel and Chen (2009), Bickel, Chen and Levina (2011), Rohe, Chatterjee and Yu (2011), Chatterjee (2012), Choi, Wolfe and Airoldi (2012), Flynn and Perry (2012), Rohe and Yu (2012), Zhao, Levina and Zhu (2012), Fishkind et al. (2013)].

We show that blockmodels for co-clustering satisfy consistency properties and remain interpretable whenever separate exchangeability holds. Exchangeability is a natural condition satisfied by many network models: it characterizes permutation invariance, implying that the ordering of nodes carries no information [Bickel and Chen (2009), Hoff (2009)]. A blockmodel is an exchangeable model in which the connection probabilities are piecewise constant. Blockmodels also provide a simplified parametric approximation in the more general nonparametric setting [Bickel, Chen and Levina (2011)].

In addition to providing oracle inequalities for blockmodel $M$-estimators corresponding to profile likelihood and least squares optimizations, we show that it is possible to identify clusterings in data—what practitioners term network communities—even when the actual generative process is far from a blockmodel. The main statistical application of our results is to enable co-clustering under model misspecification. Much effort has been devoted to the task of community detection [Newman (2006), Fortunato and Barthélemy (2007), Zhao, Levina and Zhu (2011), Fienberg (2012)], but the drawing of inferential conclusions in this setting has been limited by the need to assume a correctly specified model.

Our results imply that community detection can be understood as finding a best piecewise-constant or simple function approximation to a flexible nonparametric process. In settings where the underlying generative process is not well understood and the specification of models is thus premature, such an approach is a natural first step for exploratory data analysis. This has been likened to the use of histograms to characterize exchangeable data in nonnetwork settings [Bickel and Chen (2009)].

The article is organized as follows. In Section 2, we introduce our nonparametric setting and model. In Section 3 we present oracle inequalities for co-clustering based on blockmodel fitting. In Section 4 we give our main technical result, and discuss a concrete statistical application: quantifying how the collection of co-clusterings of the data approaches that of a generative nonparametric process. We prove our main result in Section 5, by combining a construction used to establish a theory of graph limits [Borgs et al. (2006, 2008, 2012)] with statistical learning theory results on $U$-statistics [Clémençon, Lugosi and Vayatis (2008)]. In Section 6 we illustrate our results via a simulation study, and in Section 7 we relate them to other recent work. Appendices A–C contain additional proofs and technical lemmas.

Model elicitation

Recall that fitting a blockmodel to a binary array involves partitioning it into blocks. Denote by $G=(V_{1},V_{2},E)$ a bipartite graph with edge set $E$ and vertex sets $(V_{1},V_{2})$, where assignments of vertices to $V_{1}$ or $V_{2}$ are known. For example, $V_{1}$ and $V_{2}$ might represent people and locations, with edge $(i,j)$ denoting that person $i$ frequents location $j$. See Flynn and Perry (2012) and Rohe and Yu (2012) for additional examples.

For a bipartite graph $G$ represented as a binary array $A$, the appropriate notion of exchangeability is as follows.

An array $\{A_{ij}\}_{i,j=1}^{\infty}$ of binary random variables is separately exchangeable if

for all $n=1,2,\ldots,$ all permutations $\Pi_{1},\Pi_{2}$ of for all $n=1,2,\ldots,$ all permutations $\Pi_{1},\Pi_{2}$ of $1,\ldots,n$, and all $X\in\{0,1\}^{n\times n}$.

If we identify a finite set of rows and columns of $A$ with the adjacency matrix of an observed bipartite graph $G$, then it is clear that the notion of separate exchangeability encompasses a broad class of network models. Indeed, given a single observation of an unlabeled graph, it is natural to consider the class of all models that are invariant to permutation of its adjacency matrix; see Bickel and Chen (2009) and Hoff (2009) for discussion.

The assumption of separate exchangeability is the only one we will require for our results to hold. A representation of models in this class will be given by the Aldous–Hoover theorem for separately exchangeable binary arrays.

Fix a measurable mapping $\omega\dvtx^{3}\rightarrow$. Then the following model generates an exchangeable random bipartite graph $G=(V_{1},V_{2},E)$ through its adjacency matrix $A$: {longlist}[(1)]

generate $\alpha\sim\operatorname{Uniform}(0,1)$;

The Aldous–Hoover theorem states that this representation is sufficient to describe any separately exchangeable network distribution.

Let $\{A_{ij}\}_{i,j=1}^{\infty}$ be a separately exchangeable binary array. Then there exists some $\omega\dvtx^{3}\rightarrow$, unique up to measure-preserving transformation, which generates $\{A_{ij}\}_{i,j=1}^{\infty}$.

The interpretation of the exchangeable graph model of Definition 2.2 is that each vertex has a latent parameter in $$($\xi_{i}$for vertex$i$in$V_{1}$, and$\zeta_{j}$for vertex$j$in$V_{2}$) which determines its affinity for connecting to other vertices, while$\alpha$is a network-wide connectivity parameter (nonidentifiable from a single network observation). Because$\xi$and$\zeta$are latent,$\omega(x,y)$itself is identifiable only up to measure-preserving transformation, and is hence indistinguishable from any mapping$(x,y)\mapsto\omega(\alpha,\pi_{1}(x),\pi_{2}(y))$for which$\pi_{1},\pi_{2}$are in the set$\mathcal{P}$of measure-preserving bijective maps of$$ to itself.

2 The stochastic co-blockmodel

Many popular network models can be recognized as instances of Definition 2.2. For example, Hoff, Raftery and Handcock (2002), Airoldi et al. (2008) and Kim and Leskovec (2012) all present models in which the resulting $\omega(\alpha,x,y)$ is constant in $\alpha$, while Miller, Griffiths and Jordan (2009) require the full parameterization $\omega(\alpha,x,y)$. The stochastic co-blockmodel specifies $\omega(\alpha,x,y)$ constant in $\alpha$ and also piecewise-constant in $x$ and $y$, and thus can be viewed as a simple function approximation to $\omega(x,y)$ in Definition 2.2.

Fix integers $K_{1},K_{2}>0$, a matrix $\theta\in^{K_{1}\times K_{2}}$ and discrete probability measures $\mu$ and $\nu$ on $\{1,\ldots,K_{1}\}$ and $\{1,\ldots,K_{2}\}$. Then the stochastic co-blockmodel generates an exchangeable bipartite graph $G=(V_{1},V_{2},E)$ through the matrix $A$ as follows: {longlist}[(1)]

as the mapping corresponding to Definition 2.2, with $F_{\mu}^{-1}(x)=\inf_{z}\{F_{\mu}(z)\geq x\}$ the inverse distribution function corresponding to a given distribution $\mu$.

Without loss of generality we assume $K_{1}=K_{2}=K$ in what follows, noting that our results do not depend in any crucial way on this assumption. Thus, a stochastic blockmodel’s vertices in $V_{1}$ belong to one of $K$ latent classes, as do those in $V_{2}$. Vectors $S\in\{1,\ldots,K\}^{m}$ and $T\in\{1,\ldots,K\}^{n}$ of categorical variables specify these class memberships. The matrix $\theta\in^{K\times K}$ indexes the corresponding connection affinities between classes in $V_{1}$ and $V_{2}$. Because $S$ and $T$ are latent, the stochastic co-blockmodel is identifiable only up to a permutation of its class labels.

Oracle inequalities for co-clustering

If we assume that the separately exchangeable data model of Definition 2.2 is in force, then a natural first step is to approximate $\omega(x,y)$ by way of some piecewise-constant $\omega_{\phi}(x,y)$, according to the stochastic co-blockmodel of Definition 2.3. This approximation task is equivalent to fixing $K$ and estimating $\phi=(\mu,\nu,\theta)$ by co-clustering the entries of an observed adjacency matrix $A\in\{0,1\}^{m\times n}$.

To accomplish this task, we consider $M$-estimators that involve an optimization over the latent categorical variable vectors $S\in\{1,\ldots,K\}^{m}$ and $T\in\{1,\ldots,K\}^{n}$. The resulting blockmodel estimates will reside in a set $\Phi$ containing triples $(\mu,\nu,\theta)\in\Omega_{m}\times\Omega_{n}\times^{K\times K}$, where we define $\Omega_{m}$ to be the set of all probability distributions over $\{1,\ldots,K\}$ whose elements are integer multiples of $1/m$,

and likewise for $\Omega_{n}$. Note that $\Omega_{m}$ and $\Omega_{n}$ are subsets of the standard $K-1$-simplex, chosen to contain all measures $\mu$ and $\nu$ that can be obtained by empirically co-clustering the elements of an $m\times n$-dimensional binary array. Thus, by construction, any estimator $\hat{\phi}(A)=(\hat{\mu},\hat{\nu},\hat{\theta})$ based on an empirical co-clustering of an observed binary array $A\in\{0,1\}^{m\times n}$ has codomain $\Phi$.

Given a specific $\mu$ and $\nu$, let $\mathcal{Q}_{\mu}^{m}$ denote the set of all node-to-class assignment functions that partition the set $\{1,\ldots,m\}$ into $K$ classes in a manner that respects the proportions dictated by $\mu=(\mu_{1},\ldots,\mu_{K})\in\Omega_{m}$,

and likewise for $\mathcal{Q}_{\nu}^{n}$.

2 Oracle inequalities

We now establish that, for $L^{2}$ risk and Kullback–Leibler divergence, there exist $M$-estimators that enable us to determine, with rate of convergence $n^{-1/4}$, optimal piecewise-constant approximations of the generative $\omega(x,y)$, up to quantization due to the discreteness of $\Phi$.

Let $A\in\{0,1\}^{m\times n}$ be a separately exchangeable array generated by some $\omega$ in accordance with Definition 2.2, and consider fitting a $K$-class stochastic co-blockmodel parameterized by $\phi\equiv(\mu,\nu,\theta)$ to $A$. Then as $n\rightarrow\infty$, with $K$ and $m/n$ fixed:

For the least squares co-blockmodel $M$-estimator

Given any $\phi=(\mu,\nu,\theta)$, let $B(\phi)=\max_{1\leq a,b,\leq K}|\log(\theta_{ab}/(1-\theta_{ab}))|$. Consider the profile likelihood co-blockmodel $M$-estimator

If $\phi^{*}=\operatorname{argmax}_{\phi\in\Phi}L_{\omega}(\phi)$ exists, and $B(\phi^{*})$ and $B(\hat{\phi})$ are finite, then

Theorem 3.1 can be viewed as analyzing maximum likelihood techniques in the context of model misspecification [White (1982)], and is proved in Appendix A. It establishes that minimization of the squared error between a fitted co-blockmodel and an observed binary array according to (1) serves as a proxy for approximation of $\omega$ by $\omega_{\phi}$ in mean square, and that fitting a stochastic co-blockmodel via profile likelihood according to (3.1) is equivalent to minimizing the average Kullback–Leibler divergence of the approximation $\omega_{\phi}(x,y)$ from the generative $\omega(x,y)$.

The existence of a limiting object $\omega(x,y)$ implies that we are in the dense graph regime, with expected network degree values increasing linearly as a function of $m$ or $n$. Given a correctly specified generative blockmodel, profile likelihood estimators are known to be consistent even in the sparse graph setting of polynomial or poly-logarithmic expected degree growth [Bickel and Chen (2009)]. In our setting, however, the generative model is no longer necessarily a blockmodel; in this context, both Borgs et al. (2008) and Chatterjee (2012) leave open the question of consistently estimating sparse network parameters, while Bickel, Chen and Levina (2011) give an identifiability result extending to the sparse case. The simulation study reported in Section 6 below suggests that the behavior of blockmodel estimators is qualitatively similar across at least some families of dense and sparse models.

3 Additional remarks on Theorem 3.1

In essence, Theorem 3.1 implies that the binary array $A$ yields information on its underlying generative $\omega(x,y)$ at a rate of at least $n^{-1/4}$. While the necessary optimizations in (1) and (3.1) are not currently known to admit efficient exact algorithms, they strongly resemble existing objective functions for community detection for which many authors have reported good heuristics [Newman (2006), Fortunato and Barthélemy (2007), Zhao, Levina and Zhu (2011)]. Furthermore, polynomial-time spectral algorithms are known in certain settings to find correct labelings under the assumption of a generative blockmodel [Rohe, Chatterjee and Yu (2011), Fishkind et al. (2013)], suggesting that efficient algorithms may exist when distinct clusterings or community divisions are present in the data. In this vein, Chatterjee (2012) has recently proposed a universal thresholding procedure based on the singular value decomposition.

We may replace the objective function of (3.1)with the full profile likelihood function $\max_{S\in\mathcal{Q}_{\mu}^{m},T\in\mathcal{Q}_{\nu}^{n}}\{\sum_{i=1}^{m}\log\mu_{{S(i)}}+\break\sum_{j=1}^{n}\log\nu_{{T(j)}}+\sum_{i=1}^{m}\sum_{j=1}^{n}\{A_{ij}\log\theta_{S(i)T(j)}+(1-A_{ij})\log(1-\theta_{S(i)T(j)})\}\}$. The same rate of convergence can then be established with respect to the corresponding term for $L_{\omega}(\phi)$, adapting the proofs in Appendices A and B.

Assume $\phi^{*}=\operatorname{argmax}_{\phi\in\Phi}L_{\omega}(\phi)$ exists. Terms $B(\phi^{*})$ and $B({\hat{\phi}})$ in (3) show that elements of $\theta^{*}$ and ${\hat{\theta}}$ must not approach or $1$ too quickly as $n\rightarrow\infty$; otherwise $L_{\omega}(\hat{\phi})$ can be much smaller than $L_{\omega}(\phi^{*})$.

This is a natural consequence of the fact that the Kullback–Leibler divergence of $\omega_{\phi}$ from $\omega$ is finite if and only if $\omega$ is absolutely continuous with respect to $\omega_{\phi}$. To see the implication, consider $\xi,\zeta$, and $A$ generated according to Definition 2.2 with $\omega(x,y)=1\{x\leq 1/2\}1\{y\leq 1/2\}$. Let $\mu_{1}=m^{-1}\sum_{i=1}^{m}1\{\xi_{i}\leq 1/2\}$ and $\nu_{1}=n^{-1}\sum_{j=1}^{n}1\{\zeta_{i}\leq 1/2\}$. Then the maximum-likelihood two-class blockmodel fit to $A$ will yield $\omega_{{\hat{\phi}}}(x,y)=1\{x\leq\mu_{1}\}1\{y\leq\nu_{1}\}$, and so $L_{\omega}({\hat{\phi}})$ diverges to $-\infty$ unless $\mu_{1}=\nu_{1}=1/2$.

Convergence of co-cluster estimates

We now give our main technical result and show its statistical application in enabling us to interpret the convergence of co-cluster estimates. The estimators of Theorem 3.1 require optimizations over the set of all possible co-clusterings of the data; that is, over vectors $S$ and $T$ that map the observed vertices to $1,\ldots,K$. Analogously, one may also envision an uncountable set of co-clusterings of the generative model, which map the unit interval $$to$1,\ldots,K$. We define these two sets of co-clusterings more formally and then give a result showing in what sense they become close with increasing$m$and$n$, so that optimizing over co-clusters of the data is asymptotically equivalent to optimizing over co-clusters of the generative model. This result yields the rate of convergence$\mathcal{O}_{P}(n^{-1/4})$ appearing in Theorem 3.1, and also has a geometric interpretation that sheds light on the estimators defined by (1) and (3.1).

Given a bipartite graph $G=(V_{1},V_{2},E)$ with adjacency matrix $A\in\{0,1\}^{m\times n}$, recall that the latent class vectors $S\in\{1,\ldots,K\}^{m}$ and $T\in\{1,\ldots,K\}^{n}$ respectively partition $V_{1}$ and $V_{2}$ into $K$ subsets each. To relate an empirical co-clustering of $A$ to a piecewise-constant approximation of some $\omega$, we first define the matrix $A/ST\in^{K\times K}$ to index the proportion of edges spanning each of the $K^{2}$ subset pairs defined by $S$ and $T$,

Second, we define mappings $\sigma,\tau\dvtx\rightarrow\{1,\ldots,K\}$, which will play a role analogous to $S$ and $T$. Given some $\omega\dvtx^{2}\rightarrow$, this allows us to define a matrix $\omega/\sigma\tau\in^{K\times K}$ which encodes the mass of $\omega$ assigned to each of the $K^{2}$ subset pairs defined by $\sigma$ and $\tau$ as follows:

We will use the $K\times K$ matrices $A/ST$ and $\omega/\sigma\tau$ to index all possible co-clusterings that can be induced by partitioning an observed binary array $A\in\{0,1\}^{m\times n}$ into $K^{2}$ blocks. To link these sets of co-clusters, recall from Section 3 the sets $\mathcal{Q}_{\mu}^{m}$ and $\mathcal{Q}_{\nu}^{n}$ of all node-to-class assignment functions that partition $\{1,\ldots,m\}$ and $\{1,\ldots,n\}$ into $K$ classes in manners that respect the proportions dictated by $\mu=(\mu_{1},\ldots,\mu_{K})\in\Omega_{m}$ and $\nu=(\nu_{1},\ldots,\nu_{K})\in\Omega_{n}$. Analogously, we define $\mathcal{Q}_{\mu}$ (resp., $\mathcal{Q}_{\nu}$) to be the set of partitions of $$into$K$subsets whose cardinalities are of proportions$\mu_{1},\ldots,\mu_{K}$:

We are now equipped to introduce sets ${\mathcal{F}_{\mu\nu}^{A}}$ and ${\mathcal{F}_{\mu\nu}^{\omega}}$, which describe all possible co-clusterings that can be induced from $A$ and $\omega$ with respect to $(\mu,\nu)\in\Omega_{m}\times\Omega_{n}$, and to define the related notion of a support function.

We will show below that $\sup_{\Gamma\in^{K\times K}}|h_{{{\mathcal{F}_{\mu\nu}^{A}}}}(\Gamma)-h_{{{\mathcal{F}_{\mu\nu}^{\omega}}}}(\Gamma)|$ converges in probability to zero at a rate of at least $n^{{-1/4}}$, and this result in turn gives rise to Theorem 3.1. To see why, observe that for any $(\mu,\nu,\theta)\in\Phi$, the least squares objective function of (1) can be expressed using $h_{{\mathcal{F}_{\mu\nu}^{A}}}(\theta)$ as follows:

As we prove in Appendix A, this line of argument establishes the following.

For any $(\mu,\nu,\theta)\in\Phi$, the difference between the least squares objective function of (1) and the $L^{2}$ risk $R_{\omega}$ is equal to

and the difference between the profile likelihood function of (3.1) and $L_{\omega}$ is $B(\theta)(h_{{\mathcal{F}_{\mu\nu}^{A}}}(\Gamma_{\theta})-h_{{\mathcal{F}_{\mu\nu}^{\omega}}}(\Gamma_{\theta}))$ whenever $0<\theta_{ab}<1$ for all $a,b=1,\ldots,K$, with $\Gamma_{\theta}\in^{{K\times K}}$ given element-wise by $(\Gamma_{\theta})_{ab}=\log(\theta_{ab}/(1-\theta_{ab}))/B(\theta)$.

2 A general result on consistency of co-clustering

From Lemma 4.1 we see that closeness of $h_{{{\mathcal{F}_{\mu\nu}^{A}}}}$ to $h_{{{\mathcal{F}_{\mu\nu}^{\omega}}}}$ implies closeness (up to constant terms) of the least squares objective function of (1) to the $L^{2}$ risk $R_{{\omega}}(\phi)$, and of the profile likelihood of (3.1) to the average Kullback–Leibler divergence of $\omega_{\phi}(x,y)$ from the generative $\omega(x,y)$. Equipped with this motivation, we now state our main technical result, which serves to establish the rate of convergence $\mathcal{O}_{P}(n^{{-1/4}})$ in Theorem 3.1. Its proof follows in Section 5 below.

Let $A\in\{0,1\}^{m\times n}$ be a separately exchangeable array generated by some $\omega$ in accordance with Definition 2.2. Then for each $K$ and each ratio $m/n$, there exists a universal constant $C$ such that as $n\rightarrow\infty$,

Formally, Theorem 4.1 has the following geometric interpretation:

The result of Theorem 4.1 is equivalent to the following: The Hausdorff distance between the convex hulls of ${\mathcal{F}_{\mu\nu}^{A}}$ and ${\mathcal{F}_{\mu\nu}^{\omega}}$ is $\mathcal{O}_{P}(n^{-1/4})$.

Since Theorem 4.1 holds for $\sup_{\Gamma\in^{K\times K}}|h_{{{\mathcal{F}_{\mu\nu}^{A}}}}(\Gamma)-h_{{{\mathcal{F}_{\mu\nu}^{\omega}}}}(\Gamma)|$, the leftmost inequality implies that it also holds for $\sup_{\|\Gamma\|=1}|h_{{{\mathcal{F}_{\mu\nu}^{A}}}}(\Gamma)-h_{{{\mathcal{F}_{\mu\nu}^{\omega}}}}(\Gamma)|$. Now suppose instead that Theorem 4.1 holds for $\sup_{\|\Gamma\|=1}|h_{{{\mathcal{F}_{\mu\nu}^{A}}}}(\Gamma)-h_{{{\mathcal{F}_{\mu\nu}^{\omega}}}}(\Gamma)|$; by the rightmost inequality, it then also holds for $K^{-1}\sup_{\Gamma\in^{K\times K}}|h_{{{\mathcal{F}_{\mu\nu}^{A}}}}(\Gamma)-h_{{{\mathcal{F}_{\mu\nu}^{\omega}}}}(\Gamma)|$. Thus the result of Theorem 4.1 is equivalent to the statement that

This geometric interpretation is helpful in relating our work to a series of papers by Borgs et al. (2006, 2008, 2012), which explore dense graph limits in depth and statistical applications thereof. Very broadly speaking, Borgs et al. (2008), Theorem 2.9 and Borgs et al. (2012), Theorem 4.6, analyze sets termed quotients, which resemble $\bigcup_{\mu,\nu}{\mathcal{F}_{\mu\nu}^{A}}$ and $\bigcup_{\mu,\nu}{\mathcal{F}_{\mu\nu}^{\omega}}$. The authors show convergence of these sets in the Hausdorff metric at rate $\mathcal{O}(\log^{{-1/2}}n)$, based on a distance termed the cut metric, and detail implications that can also be related to those of Bickel, Chen and Levina (2011).

In fixing $\mu$ and $\nu$ through our $M$-estimators, we are studying what Borgs et al. term the microcanonical quotients. Because our results require only convergence of the closed convex hulls of ${\mathcal{F}_{\mu\nu}^{A}}$ and ${\mathcal{F}_{\mu\nu}^{\omega}}$, we are able to obtain an exponentially faster bound on the rate of convergence.

3 Interpreting convergence of blockmodel estimates

Recall that the $M$-estimators of Theorem 3.1 each involve an optimization over the set ${\mathcal{F}_{\mu\nu}^{A}}$ by way of its support function, which in turn represents its convex hull. Suppose that ${\hat{\phi}}\equiv(\hat{\mu},\hat{\nu},{\hat{\theta}})$ optimizes either objective function in Theorem 3.1. Then the following corollary of Theorem 3.1 shows that $\hat{\phi}$ is interpretable, in that there will exist a partition $\hat{\sigma},\hat{\tau}$ of $\omega$ yielding co-clusters of equal size and asymptotically equivalent connectivity.

Let $\hat{\phi}=(\hat{\mu},\hat{\nu},\hat{\theta})$ minimize the least squares criterion of (1). Then there exists some pair $(\hat{\sigma},\hat{\tau})\in\mathcal{Q}_{\hat{\mu}}\times\mathcal{Q}_{\hat{\nu}}$ such that

Similarly, if $\hat{\phi}=(\hat{\mu},\hat{\nu},\hat{\theta})$ maximizes the profile likelihood criterion of (3.1) and $\phi^{*}=\operatorname{argmax}_{\phi\in\Phi}L_{\omega}(\phi)$ exists, then there is some $(\hat{\sigma},\hat{\tau})\in\mathcal{Q}_{\hat{\mu}}\times\mathcal{Q}_{\hat{\nu}}$ with

where $D(p\|p^{\prime})=p\log(p/p^{\prime})+(1-p)\log[(1-p)/(1-p^{\prime})]\geq 0$ is the Kullback–Leibler divergence of a $\operatorname{Bernoulli}(p^{\prime})$ distribution from a $\operatorname{Bernoulli}(p)$ one.

We show the latter result; parallel arguments yield the former. Since $\omega_{\hat{\phi}}(x,y)={\hat{\theta}}_{{F^{-1}_{\hat{\mu}}(x)F^{-1}_{\hat{\nu}}(y)}}$ for the co-blockmodel, by letting $\sigma$ and $\tau$ satisfy $\sigma(x)=F^{{-1}}_{\hat{\mu}}(\pi_{1}(x))$ and $\tau(y)=F^{{-1}}_{\hat{\nu}}(\pi_{2}(y))$ we may express $L_{\omega}({\hat{\phi}})$ as

Thus, for any $\varepsilon>0$, there exists some choice of $(\hat{\sigma},\hat{\tau})\in\mathcal{Q}_{{\hat{\mu}}}\times\mathcal{Q}_{{\hat{\nu}}}$ such that

If we now take $\hat{\theta}^{(\omega)}_{{ab}}=(\omega/\hat{\sigma}\hat{\tau})_{ab}/(\hat{\mu}_{a}\hat{\nu}_{b})$ for $a,b=1,\ldots,K$, we see by a similar argument that since $L_{\omega}(\phi^{*})=\max_{\phi\in\Phi}L_{\omega}(\phi)$, we have in turn that

Expanding $D(\hat{\theta}^{(\omega)}_{ab}\|\hat{\theta}_{ab})$ in accordance with its definition, we then see that

Choosing $\varepsilon=o(n^{{-1/4}})$ and applying Theorem 3.1 completes the proof.

Corollary 4.2 ensures that co-blockmodel fits remain interpretable, even in the setting of model misspecification. It establishes that the identification of co-clusters in an observed exchangeable binary array $A$ indicates with high probability the existence of co-clusters of equal size and asymptotically equivalent connectivity in the underlying generative process $\omega$.

Proof of Theorem 4.1

Our proof strategy is inspired by Borgs et al. (2008) and adapts certain of its tools, but also requires new techniques in order to attain polynomial rates of convergence. Most significantly, we do not use the Szemerédi regularity lemma, which typically features strongly in the graph-theoretic literature, and provides a means of partitioning any large dense graph into a small number of regular clusters. Results in this direction are possible, but instead we use a Rademacher complexity bound for $U$-statistics adapted from Clémençon, Lugosi and Vayatis (2008), allowing us to achieve the improved rates of convergence described above.

The main step in proving Theorem 4.1 is to establish pointwise convergence of $h_{{\mathcal{F}_{\mu\nu}^{A}}}(\Gamma)$ to $h_{{\mathcal{F}_{\mu\nu}^{\omega}}}(\Gamma)$ for any fixed $\Gamma$. We do this through Proposition 5.1 below, after which we may apply it to a union bound over a covering of all $\Gamma\in^{K\times K}$ to deduce the result of Theorem 4.1. Appendix B provides a formal statement and proof of this argument, along with proofs of all supporting lemmas.

Assume the setting of Theorem 4.1, fixing $m=\rho n$. Then there exist constants $C_{K},n_{K}$ such that, given any $\Gamma\in^{K\times K}$, $\mu$, $\nu$, $\omega$, and $A\in\{0,1\}^{m\times n}$ generated from $\omega$, it holds for all $n\geq n_{K}$ that

To obtain the claimed result, we must establish lower and upper bounds on the support function $h_{{\mathcal{F}_{\mu\nu}^{A}}}(\Gamma)$ that show its convergence to $h_{{\mathcal{F}_{\mu\nu}^{\omega}}}(\Gamma)$ at rate $\mathcal{O}_{P}(n^{{-1/4}})$. Recalling the definitions of $h_{{\mathcal{F}_{\mu\nu}^{A}}}(\Gamma)$ and $h_{{\mathcal{F}_{\mu\nu}^{\omega}}}(\Gamma)$ in (4), we first require a statement of Lipschitz conditions on $\langle\Gamma,A/ST\rangle$ and $\langle\Gamma,\omega/\sigma\tau\rangle$. Its proof follows by direct inspection.

Define for measurable mappings $\sigma,\sigma^{\prime}$ over $$ the metric

and analogously the standard Hamming distance for sequences, with respect to normalized counting measure. Then for any $\Gamma\in^{K\times K}$ and $A,A^{\prime}\in^{m\times n}$, with $(S,T,\omega,\sigma,\tau)$ as defined in Section 4.1, we have that: {longlist}[(1)]

$|\langle\Gamma,A/ST\rangle-\langle\Gamma,A^{\prime}/ST\rangle|\leq 1/(mn)$ if $A,A^{\prime}$ differ by a single entry.

In conjunction with McDiarmid’s inequality, these Lipschitz conditions yield the following lower bound on $h_{{\mathcal{F}_{\mu\nu}^{A}}}(\Gamma)$, proved in Appendix B.1.

Assume the setting of Theorem 4.1. Then there exist constants $C{}^{\prime}_{K},n{}^{\prime}_{K}$ such that, given any $\Gamma\in^{K\times K},\mu,\nu,\omega$, and $A\in\{0,1\}^{\rho n\times n}$ generated from $\omega$, for all $n\geq n_{K}^{\prime}$,

The upper bound comes by way of Rademacher complexity arguments. The remainder of this section and Appendix B is devoted to its proof.

Assume the setting of Theorem 4.1. Then there exist constants $C^{\prime\prime}_{K},n^{\prime\prime}_{K}$ such that, given any $\Gamma\in\break^{K\times K},\mu,\nu,\omega$ and $A\in\{0,1\}^{\rho n\times n}$ generated from $\omega$, for all $n\geq n_{K}^{\prime\prime}$,

Proposition 5.1 now follows simply by combining Lemmas 5.2 and 5.3.

Lemma 5.3 represents the main technical hurdle in obtaining the polynomial rate of convergence given in Theorems 3.1 and 4.1. To illustrate the main ideas as clearly as possible, we will introduce our Rademacher complexity arguments below for the case $K=2$, deferring the necessary generalizations to Appendix B.

We first define $W\in^{m\times n}$ with reference to Definition 2.2 as

and then define, in direct analogy to $h_{{\mathcal{F}_{\mu\nu}^{A}}}(\Gamma)$,

Fix some measurable $\omega\dvtx^{2}\rightarrow$, with $W\in^{m\times n}$ generated by $\omega$ and $A\in\{0,1\}^{m\times n}$ generated by $W$, and some $\Gamma\in^{K\times K}$. Then for any $\varepsilon>0$,

Let $\mathcal{I}$ and $\mathcal{J}$ be sets of deterministic size, whose elements are sampled without replacement from $1,\ldots,m$ and $1,\ldots,n$. Let $W$ be generated as in Lemma 5.4, and fix $\Gamma\in^{K\times K}$. Given $W,\mathcal{I},\mathcal{J}$ and $(Q,R)\in\mathcal{Q}_{\mu}^{m}\times\mathcal{Q}_{\nu}^{n}$, let $\hat{S}^{R}\equiv\hat{S}^{R,\mathcal{J},W}$ and $\hat{T}^{Q}\equiv\hat{T}^{Q,\mathcal{I},W}$ denote partitions satisfying

To bound the right-hand side of (5.5) relative to $h_{{\mathcal{F}_{\mu\nu}^{\omega}}}(\Gamma)$, we will introduce an additional construction comprising several steps. Specifically, for fixed $(Q,R)$ and $\Gamma$, we will define function classes $\mathcal{Q}_{U}$ and $\mathcal{Q}_{V}$, and a random functional $G_{\sigma\tau}$ which approximates $\langle\Gamma,W/\hat{S}^{R}\hat{T}^{Q}\rangle$ for some $(\hat{\sigma},\hat{\tau})\in\mathcal{Q}_{U}\times\mathcal{Q}_{V}$. By a Rademacher complexity argument, $G_{\hat{\sigma}\hat{\tau}}$ will concentrate for all $(Q,R)$ near its expectation, which itself will be bounded by $h_{\mathcal{F}_{\mu\nu}^{\omega}}(\Gamma)$.

and so $\hat{S}^{R}$ will assign to class $1$ the $\mu_{1}m$ largest elements of $U(\xi_{1}),\ldots,U(\xi_{m})$. If $U$ is invertible, this set can be written $\{\xi_{i}\dvtx U(\xi_{i})<t\}$ for some $t$. To treat noninvertible $U$, define $\mathcal{Q}_{U}$ to be the class of functions $\{1_{u}\dvtx u\in\}$, with $1_{u}$ a one-sided interval on the range of $U$ with lexicographic “tie-breaking”:

Then there exists $\hat{\sigma}\in\mathcal{Q}_{U}$ such that $\hat{S}^{R}$ can be chosen to satisfy

Let $V$ denote a function defined analogously to $U$ as follows:

and likewise define $\mathcal{Q}_{V}$ so that there exists $\hat{\tau}\in\mathcal{Q}_{V}$ such that $\hat{T}^{Q}$ can be chosen to satisfy

We are now ready to define $G_{\sigma\tau}$. Given any $\sigma\in\mathcal{Q}_{U}$ and $\tau\in\mathcal{Q}_{V}$, let

where $\overline{\mathcal{I}}$ is the complement of $\mathcal{I}$ in $\{1,\ldots,m\}$, and $\overline{\mathcal{J}}$ the complement of $\mathcal{J}$ in $\{1,\ldots,n\}$. Comparing $G_{\sigma\tau}$ to Lemma 5.5, we see that $G_{\hat{\sigma}\hat{\tau}}$ well approximates $\langle\Gamma,W/\hat{S}^{R}\hat{T}^{Q}\rangle$ whenever $|\mathcal{I}|$ and $|\mathcal{J}|$ are small; and indeed, we will later set $|\mathcal{I}|=|\mathcal{J}|=n^{1/2}$ in order to obtain an upper bound for $h_{{\mathcal{F}_{\mu\nu}^{A}}}(\Gamma)-h_{\mathcal{F}_{\mu\nu}^{\omega}}(\Gamma)$.

By construction, the random classes $\mathcal{Q}_{U}$ and $\mathcal{Q}_{V}$ are independent of the random variables $\{\xi_{i}\}_{i\in\overline{\mathcal{I}}}$ and $\{\zeta_{j}\}_{i\in\overline{\mathcal{J}}}$ appearing in the summand of $G_{\sigma\tau}$. As a result, we may bound the deviation $\delta_{UV}$ of $G_{\sigma\tau}$ from its expectation,

using Rademacher complexity results for $U$-statistics due to Hoeffding (1963) and Clémençon, Lugosi and Vayatis (2008), Lemma A.1, applied to the class of one-sided interval functions.

Lemma 5.6 is proved in Appendix B.5 to hold for arbitrary $K$, under the appropriate generalization of $\mathcal{Q}_{U},\mathcal{Q}_{V}$, and quantities that depend on them.

Similarly, we may bound $\delta_{U}$, defined for $K=2$ as the maximum discrepancy between the expected and empirical class frequency in $\mathcal{Q}_{U}$,

with $\delta_{V}$ defined mutatis mutandis. We then have the following result, proved for arbitrary $K$ (with appropriate redefinitions of $\delta_{U},\delta_{V}$) in Appendix B.6.

We state and prove a final auxiliary lemma prior to the proof of Lemma 5.3.

after which we may upper-bound terms (i)–(iv) in turn as follows.

First, since $|\omega(x,y)\Gamma_{\hat{\sigma}(x)\hat{\tau}(y)}|\leq 1$ for all $(x,y)$, it follows from their respective definitions that $\langle\Gamma,W/\hat{S}^{R}\hat{T}^{Q}\rangle-G_{\hat{\sigma}\hat{\tau}}(\xi,\zeta)$ is deterministically bounded above by $|\mathcal{I}|/m+|\mathcal{J}|/n$. Hence, term (i) is bounded by the same quantity.

To conclude, note term (iv) is deterministically upper-bounded by .

We may now establish the claimed upper bound on $h_{{\mathcal{F}_{\mu\nu}^{A}}}(\Gamma)-h_{\mathcal{F}_{\mu\nu}^{\omega}}(\Gamma)$.

Proof of Lemma 5.3 Combining the results of Lemmas 5.4–5.8 yields directly that, with probability at least $1-2e^{-2mn\varepsilon^{2}/(m+n)}-2K^{m+n}e^{-2mn\varepsilon^{2}}$,

with probability at least $1-2e^{-\sqrt{n}[2\rho/(\rho+1)]}-2K^{(\rho+1)n}e^{-2\rho n^{3/2}}$. Thus we have established the claimed upper bound on $h_{\mathcal{F}_{\mu\nu}^{A}}(\Gamma)$ in terms of $h_{\mathcal{F}_{\mu\nu}^{\omega}}(\Gamma)$.

Simulation study

We now present a brief simulation study which investigates empirical rates of convergence as model misspecification increases. We control the degree of misspecification through a sigmoidal functional form $f_{\beta}(x)\dvtx\rightarrow[-1/2,1/2]$, parameterized by $\beta\geq 1$,

Each $f_{\beta}(x)$ describes a strictly monotone increasing sigmoidal curve on $$, proportional to$x-1/2$for$\beta=1$and to$1\{x>1/2\}-1/2$in the limit as$\beta\rightarrow\infty$. Normalization by$Z_{\beta}$maintains constant area under$|f_{\beta}|$.

To explore sparse graph regimes, we introduce an additional $n$-dependent parameter $\rho_{n}\in(0,1)$, and take the outer product $f_{\beta}(x)f_{\beta}(y)$ to obtain a separable generative function $\rho_{n}\omega_{\beta}(x,y)=\rho_{n}(f_{\beta}(x)f_{\beta}(y)+1/2)$. As $\beta\rightarrow\infty$, this tends to a stochastic co-blockmodel, with two classes of equal size.

Figure 1 shows a number of simulation results based on this model. Specifically, for $\beta\in\{1,3,5\}$ and $\rho_{n}\in\{0.5,n^{-2/3},n^{-1}\log^{2}n\}$, one thousand separable $n\times n$ binary arrays were generated from the corresponding $\rho_{n}\omega_{\beta}(x,y)$, for network sizes ranging from 100–500 for dense graphs (left column), and 100–2200 for sparse graphs (right columns). We see immediately that the simulation results of Figure 1 are qualitatively similar for all three regimes, suggesting that at least in some cases, co-blockmodel estimators will converge despite model misspecification in sparse as well as dense graph regimes.

Each of the $n\times n$ arrays described above was fitted by a two-class co-blockmodel, whose parameters $\hat{\phi}=(\hat{\mu},\hat{\nu},\hat{\theta})$ were obtained by heuristically optimizing the profile likelihood criterion of (3.1) using an algorithmic approach based on simulated annealing [Choi, Wolfe and Airoldi (2012)]. Parameter values were initialized to coincide with the optimal blockmodel approximation based on $\rho_{n}\omega_{\beta}/\sigma^{*}\tau^{*}$, where $\sigma^{*},\tau^{*}\dvtx\rightarrow\{1,2\}$ each map the interval $[0,1/2)$ to class $1$ and the interval $[1/2,1]$ to class $2$.

Lemma C.1 establishes that $\phi^{*}=\operatorname{argmax}_{\phi\in\Phi}L_{\rho_{n}\omega_{\beta}}(\phi)$ exists in this setting, and that $L_{\rho_{n}\omega_{\beta}}(\phi)$ may be straightforwardly computed for any triple $\phi=(\mu,\nu,\theta)$ of two-class co-blockmodel parameters. Corollary C.1 then yields a finite set containing $\phi^{*}=\operatorname{argmax}_{\phi\in\Phi}L_{\rho_{n}\omega_{\beta}}(\phi)$, from which we found that $\phi^{*}$ corresponded to the blockmodel induced by $\sigma^{*}$ and $\tau^{*}$. Thus we were able to evaluate the relative excess risk $[L_{\rho_{n}\omega_{\beta}}(\phi^{*})-L_{\rho_{n}\omega_{\beta}}(\hat{\phi})]/L_{\rho_{n}\omega_{\beta}}(\phi^{*})$, shown as a percentage in the top row of Figure 1, and seen to decay toward .

The bottom row of Figure 1 shows the normalized Kullback–Leibler divergences $\rho_{n}^{-1}D(\rho_{n}\omega_{\beta}\|\rho_{n}\omega_{{\hat{\phi}}})$ decaying toward the grey horizontal lines representing the limiting values of $D(\rho_{n}\omega_{\beta}\|\rho_{n}\omega_{{\phi^{*}}})$ as $\rho_{n}\rightarrow 0$. These are order-one quantities, obtained through a Taylor expansion of $D(\rho_{n}\omega_{\beta}\|\rho_{n}\omega_{{\phi^{*}}})$. Smaller divergences are achieved when $\beta$ is large, reflecting the fact that as $\beta$ increases, $\rho_{n}\omega_{\beta}(x,y)$ becomes closer to a co-blockmodel.

Overall, we see that the simulation results shown in Figure 1 are consistent with the behavior predicted by Theorem 3.1 for profile likelihood maximization; qualitatively similar results were also obtained for the least squares setting of Theorem 3.1 and hence are omitted for brevity.

Discussion

In this article we have addressed the case of network co-clustering, in which the inference task is to group two sets of network nodes into classes based on their observed relations. Our results significantly generalize known consistency results for the blockmodel and its co-blockmodel variant: they do not require the data to be generated (even approximately) by a co-blockmodel, and they achieve improved rates of convergence relative to results from the graph limits literature, through the use a Rademacher complexity bound for $U$-statistics adapted from Clémençon, Lugosi and Vayatis (2008). The assumption of a nonparametric generative model is both more general and more realistic, and to our knowledge Theorems 3.1 and 4.1 are the first for this regime to establish polynomial rates of convergence.

In the work of Clémençon, Lugosi and Vayatis (2008), these Rademacher complexity results are used to derive convergence rates for learning pairwise rankings. This setting is related to ours, but differs in some important ways. Those authors seek a rule $r\dvtx\mathcal{X}\times\mathcal{X}\rightarrow\{-1,+1\}$ such that, given $X,X^{\prime}\in\mathcal{X}$, $r$ indicates which has the higher rank. In this setting, $X$ and $X^{\prime}$ can be thought of as covariates describing the two objects for which a relative ranking is desired, and $\mathcal{X}$ represents the space of allowable covariate values. In our network setting, the nonparametric model $\omega\dvtx^{2}\rightarrow$ is analogous to a ranking rule, with $\mathcal{X}$ taken to be $$. However,$X$and$X^{\prime}$ are never observed in the data, and effectively must be imputed up to measure-preserving transformation.

The recent work of Flynn and Perry (2012) analyzes the consistency of co-clustering with model misspecification, but in a rather different setting, with the data matrix $A$ assumed to be real valued, along with a real-valued generalization of the co-blockmodel. This generalization utilizes discrete latent class variables $S$ and $T$; conditioned on $S(i)$ and $T(j)$, the distribution of $A_{ij}$ is assumed to have mean $\theta_{S(i)T(j)}$, but may otherwise be arbitrary up to technical conditions, and may be misspecified in the estimator. Under these assumptions, it is shown that the latent classes can be estimated consistently if their number is known. In the case where $A$ is binary, the conditions of Flynn and Perry (2012) are equivalent to assuming a generative co-blockmodel with a known number of classes.

Finally, the very recent work of Chatterjee (2012) derives a simple and elegant spectral method to consistently estimate the matrix $W$ defined in the proof of Lemma 5.3 in Section 5.2, that is, the mapping $\omega(x,y)$, evaluated at the values of the latent variables $\xi_{1},\ldots,\xi_{m}$, and $\zeta_{1},\ldots,\zeta_{n}$. This implies consistency of estimation of $\omega$ in the $L^{2}$ sense, and while rates of convergence are not given for general $\omega$, they can be established for particular instances, such as under the assumption of a generative blockmodel whose number of classes $K$ is growing with $n$. Our setting is distinct, in that we desire only the best blockmodel approximation to $\omega$, and so are able to establish $L^{2}$ rates of convergence that are independent of $\omega$.

Appendix A Proof of Theorem 3.1 and Lemma 4.1

To prove Theorem 3.1, we first denote the objective functions of (1) and (3.1) by $R_{A}(\phi)$ and $L_{A}(\phi)$, respectively. Lemma 4.1, proved below, relates $R_{A}(\phi)-R_{\omega}(\phi)$ and $L_{A}(\phi)-L_{\omega}(\phi)$ to the support functions $h_{{\mathcal{F}_{\mu\nu}^{A}}}(\cdot)$ and $h_{{\mathcal{F}_{\mu\nu}^{\omega}}}(\cdot)$, after which the result follows directly from Theorem 4.1.

To see this, let $\hat{\phi}\equiv(\hat{\mu},\hat{\nu},\hat{\theta})=\operatorname{argmin}_{\phi\in\Phi}R_{A}(\phi)$. For any $\phi\in\Phi$, we have

where the first inequality holds because $R_{A}(\hat{\phi})-R_{A}(\phi)\leq 0$, and the second holds by the triangle inequality and Lemma 4.1. Applying Theorem 4.1 and choosing $\phi$ to satisfy $R_{\omega}(\phi)\leq\inf_{\phi^{\prime}\in\Phi}R_{\omega}(\phi^{\prime})+n^{-1/4}$ then yields the result.

Now, assume $\phi^{*}=\operatorname{argmax}_{\phi\in\Phi}L_{\omega}(\phi)$ exists, and set $\hat{\phi}=\operatorname{argmax}_{\phi\in\Phi}L_{A}(\phi)$. Whenever $0<{\hat{\theta}}_{ab},{\theta^{*}}_{ab}<1$ for all $a,b=1,\ldots,K$, the second result $[L_{\omega}(\phi^{*})-L_{\omega}(\hat{\phi})]/[B(\theta^{*})+B(\hat{\theta})]=\mathcal{O}_{P}(n^{-1/4})$ of Theorem 3.1 follows similarly from

Proof of Lemma 4.1 We show the results of the lemma directly,

where the second line follows from the definition of ${\mathcal{F}_{\mu\nu}^{A}}$, and the last line from that of $h_{{\mathcal{F}_{\mu\nu}^{A}}}$. Letting $(\sigma,\tau)$ satisfy $\sigma(x)=F^{{-1}}_{{\mu(\pi_{1}(x))}}$ and $\tau(y)=F^{{-1}}_{{\nu(\pi_{2}(y))}}$,

Following similar steps, we show the second result as follows:

since $\max_{F\in{\mathcal{F}_{\mu\nu}^{A}}}\sum_{a,b}F_{ab}B(\theta)(\Gamma_{\theta})_{ab}=B(\theta)h_{{\mathcal{F}_{\mu\nu}^{A}}}(\Gamma_{\theta})$, and similarly

Appendix B Auxiliary proofs for Theorem 4.1

Below we provide proofs of all supporting lemmas for Theorem 4.1, and state and prove the covering argument used to establish the theorem: {longlist}[(1)]

First, in Sections B.1–B.3 below, we prove auxiliary Lemmas 5.2, 5.4 and 5.5 as stated in Section 5.

Then, in Section B.4, we generalize the definitions of $\mathcal{Q}_{U}$ and $\mathcal{Q}_{V}$, given in Section 5.2 for $K=2$, to arbitrary $K$; this induces generalizations of the quantities $\delta_{U},\delta_{V}$ and $\delta_{UV}$ in the natural way.

Then, in Sections B.5 and B.6, we prove Lemmas 5.6 and 5.7, which depend on $(\mathcal{Q}_{U},\mathcal{Q}_{V},\delta_{U},\delta_{V},\delta_{UV})$ as defined for arbitrary $K$.

Finally, in Section B.7, we extend the pointwise convergence result of Proposition 5.1 by way of a covering argument for all $\Gamma\in^{K\times K}$.