# Consistency and Asymptotic Normality of Latent Blocks Model Estimators

Vincent Brault, Christine Keribin, Mahendra Mariadassou

► To cite this version:

Vincent Brault, Christine Keribin, Mahendra Mariadassou. Consistency and Asymptotic Normality of Latent Blocks Model Estimators. Electronic Journal of Statistics, 2020, 14 (1), pp.1234-1268. 10.1214/20-EJS1695. hal-01511960v3

HAL Id: hal-01511960

https://hal.science/hal-01511960v3

Submitted on 25 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Consistency and Asymptotic Normality of Latent Block Model Estimators

Vincent Brault

Univ. Grenoble Alpes, CNRS, Grenoble INP*, LJK, 38000 Grenoble, France

e-mail: vincent.brault@univ-grenoble-alpes.fr

url: https://www-ljk.imag.fr/membres/Vincent.Brault/

and

Christine Keribin

Université Paris-Saclay, CNRS, Inria, Laboratoire de mathématiques d'Orsay, 91405,

Orsay, France

e-mail: christine.keribin@math.u-psud.fr

url: https://www.math.u-psud.fr/~keribin/

and

Mahendra Mariadassou

MaIAGE, INRAE, Université Paris-Saclay, 78352 Jouy-en-Josas, France

e-mail: mahendra.mariadassou@inrae.fr

url: https://mahendra-mariadassou.github.io/

**Abstract:** The Latent Block Model (LBM) is a model-based method to cluster simultaneously the *d* columns and *n* rows of a data matrix. Parameter estimation in LBM is a difficult and multifaceted problem. Although various estimation strategies have been proposed and are now well understood empirically, theoretical guarantees about their asymptotic behavior is rather sparse and most results are limited to the binary setting. We prove here theoretical guarantees in the valued settings. We show that under some mild conditions on the parameter space, and in an asymptotic regime where log(*d*)/*n* and log(*n*)/*d* tend to 0 when *n* and *d* tend to infinity, (1) the maximum-likelihood estimate of the complete model (with known labels) is consistent and (2) the log-likelihood ratios are equivalent under the complete and observed (with unknown labels) models. This equivalence allows us to transfer the asymptotic consistency, and under mild conditions, asymptotic normality, to the maximum likelihood estimate under the observed model. Moreover, the variational estimator is also consistent and, under the same conditions, asymptotically normal.

**Keywords and phrases:** Latent Block Model, asymptotic normality, Maximum Likelihood Estimate, Concentration Inequality.

## 1. Introduction

Co-clustering is an unsupervised method to cluster simultaneously the *n* rows and *d* columns of a rectangular data matrix. The assignments of each row to

<footnote>*Institute of Engineering Univ. Grenoble Alpes</footnote>
one of the row-clusters and of each column to one of the column-clusters are unknown and the aim is to determine them. Then, rows and columns can be re-ordered according to their assignments, highlighting the natural structure of the data with distinct blocks having homogeneous observations. This leads to a parsimonious data representation, as can be shown on Figure 1.

<figure><img src="image_2.png" /><figcaption>FIG 1. A binary data matrix before (left) and after (middle) row and column reordering, and its parsimonious data representation (right).</figcaption></figure>

Co-clustering can be used in numerous applications, and especially ones with large data sets, such as recommendation systems (to discover a segmentation of customers with regard to a segmentation of products), genomics (to simultaneously define groups of genes having the same expression with regards to groups of experimental conditions) or text mining (to define simultaneously groups of texts and groups of words).

Among the co-clustering methods, the Latent Block Model (LBM) defines a probabilistic model as a mixture model with latent rows and columns assignments. LBM can deal with binary [6], Gaussian [8], categorical [9] or count [7] data. Due to the complex dependency structure induced by this modeling, neither the likelihood, nor the distribution of the assignments conditionally to the observations needed in the E-step of the EM algorithm, traditionally used for mixture models, are numerically tractable. Estimation can be however performed either with a variational approximation leading to an approximate value of the maximum likelihood estimator, or with a Bayesian approach (VBayes algorithm or Gibbs sampler). For example, [9] recommends using a Gibbs sampler combined with a VBayes algorithm.

The asymptotics of the maximum likelihood (MLE) and variational (VE) estimators also raise interesting theoretical questions. This topic was first addressed for the stochastic block model (SBM) [13], where the data is a random graph encoded by its adjacency binary matrix: the rows and columns represent the nodes, so that there is only one partition, shared by rows and columns, and a unique asymptotic direction.

For a binary SBM and under the true parameter value, Theorem 3 of [4] states that the distribution of the assignments conditionally to the observations converges to a Dirac of the real assignments. Moreover, this convergence remains valid under the estimated parameter value, assuming that this estimator converges at rate at least $n^{-1}$, where n is the number of nodes (Proposition 3.8). This assumption is not trivial, and it was not established that such an estimator exists except in some particular cases ([1] for example). [10] presented a unified
frame for LBM and SBM in case of valued observations satisfying a concentration inequality, and showed the consistency of the conditional distribution of the assignments under all parameter values in a neighborhood of the true value. [3] and [2] proved the consistency and asymptotic normality of the MLE for the binary SBM but failed to account for complications induced by symmetries in the parameter. Building upon the work from [4], they first studied the asymptotic behavior of the MLE in the complete model (observations and assignments) with binary observations which is simple to handle; then, they showed that the complete likelihood and the marginal likelihood have similar asymptotic behaviors by the use of a Bernstein inequality for bounded observations.

Following the main ideas of [2], we prove that the observed likelihood ratio and the complete likelihood ratio computed at the true assignments are asymptotically equivalent, up to a multiplicative term. This term depends on some model symmetry and was omitted in [2] although it is necessary to prove the asymptotic results. We then settle the asymptotic normality of the maximum likelihood and variational estimators. All these results are stated not only for binary observations, but also more generally for observations coming from univariate exponential families in canonical form, which is essential regarding the LBM usages. This leads us to develop a Bernstein-type inequality for sub-exponential variables as the Hoeffing's concentration inequality used in [2] is only relevant for upper-bounded observations.

The paper is organized as follows. The model, main assumptions and notations are introduced in Section 2, where the concept of model symmetry is also discussed. Section 3 proves the asymptotic normality of the complete likelihood estimator, and section 4 studies conditional and profile log-likelihoods. Our main result showing that the observed likelihood ratio behaves like the complete likelihood ratio is stated in section 5, and its consequences in terms of consistency and asymptotic normality of the MLE and variational estimators are presented in section 6. Most of the proofs are postponed to the appendices to improve the general readability: appendix A for properties of conditional and profile log-likelihoods, B for the steps of the main result, C for concentration inequalities for specific sub-exponential variables and D for other technical results.

## 2. Model, assumptions and definitions

We observe a data matrix $X = (x_{ij})$ with $n$ rows and $d$ columns. The LBM assumes that there exists a latent structure in the form of the Cartesian product of a partition of $g$ row-clusters by a partition of $m$ column-clusters with the following characteristics:

* the latent row assignments $z = (z_1, \dots, z_n)$ are independent and identically distributed with a common multinomial distribution on $g$ categories:

$$z_i \sim \mathcal{M}(1, \pi = (\pi_1, \dots, \pi_g))$$

For $k = 1, \dots, g$, $z_{ik} = 1$ if row $i$ belongs to row-group $k$, 0 otherwise.
In the same way, the latent column assignments $\mathbf{w} = (\mathbf{w}_1, \dots, \mathbf{w}_d)$ are i.i.d. multinomial variables with $m$ categories:

$$ \mathbf{w}_j \sim \mathcal{M}(1, \rho = (\rho_1, \dots, \rho_m)) $$

For $\ell = 1, \dots, m$, $w_{j\ell} = 1$ if column $j$ belongs to column-group $\ell$ and 0 otherwise.

* the row and column assignments are independent: $p(\mathbf{z}, \mathbf{w}) = p(\mathbf{z})p(\mathbf{w})$

* conditionally to row and column assignments $\mathbf{z} \times \mathbf{w}$, the observed data $X_{ij}$ are independent, and their conditional distribution $\varphi(\cdot, \alpha)$ belongs to the same parametric family, which parameter $\alpha$ only depends on the given block:

$$ X_{ij} | \{z_{ik} w_{j\ell} = 1\} \sim \varphi(\cdot, \alpha_{k\ell}). $$

Hence, the complete parameter set is $\theta = (\pi, \rho, \alpha) \in \Theta$, with $\alpha = (\alpha_{11}, \dots, \alpha_{gm})$ and $\Theta$ the parameter space. Figure 2 summarizes these notations.

*Remark 2.1.* Group, class and cluster in one hand, label and assignment in the other hand will be used indistinctly. Moreover, for notation convenience, $\sum_i$, $\sum_j$, $\sum_k$, $\sum_\ell$ stand for $\sum_{i=1}^n$, $\sum_{j=1}^d$, $\sum_{k=1}^g$, $\sum_{\ell=1}^m$.

<figure><img src="image_3.png" /><figcaption>FIG 2. Notations. Left: Notations for the elements of observed data matrix are in black, notations for the block clusters are in blue. Right: Notations for the model parameter.</figcaption></figure>

When performing inference from data, we denote $\theta^* = (\pi^*, \rho^*, \alpha^*)$ the true parameter set, i.e. the parameter values used to generate the data, and $\mathbf{z}^*$ and $\mathbf{w}^*$ the true (and usually unobserved) row and column assignments. For indicator membership variables $\mathbf{z}$ and $\mathbf{w}$, we also denote:

* $z_{+k} = \sum_i z_{ik}$ and $w_{+\ell} = \sum_j w_{j\ell}$
• $z_{+k}^*$ and $w_{+\ell}^*$ their counterpart for $z^*$ and $w^*$.

The confusion matrix allows one to compare the partitions.

**Definition 2.2** (confusion matrices). For given assignments $z$ and $z^*$ (resp. $w$ and $w^*$), we define the confusion matrix between $z$ and $z^*$ (resp. $w$ and $w^*$), denoted $\mathbb{R}_g(z)$ (resp. $\mathbb{R}_m(w)$), as follows:

$$ \mathbb{R}_g(z)_{kk'} = \frac{1}{n} \sum_i z_{ik}^* z_{ik'} \quad \text{and} \quad \mathbb{R}_m(w)_{\ell\ell'} = \frac{1}{d} \sum_j w_{j\ell}^* w_{j\ell'} $$

## 2.1. Likelihood

When the labels are known, the complete log-likelihood is given by:

$$ \begin{aligned} \mathcal{L}_c(z, w; \theta) &= \log p(x, z, w; \theta) \\ &= \log \left\{ \left( \prod_{i,k} \pi_k^{z_{ik}} \right) \left( \prod_{j,\ell} \rho_{\ell}^{w_{j\ell}} \right) \left( \prod_{i,j,k,\ell} \varphi(x_{ij}; \alpha_{k\ell})^{z_{ik}w_{j\ell}} \right) \right\} \quad (2.1) \\ &= \log \left\{ \left( \prod_i \pi_{z_i} \right) \left( \prod_j \rho_{w_j} \right) \left( \prod_{i,j} \varphi(x_{ij}; \alpha_{z_i w_j}) \right) \right\}. \end{aligned} $$

In an unsupervised setting, the labels are unobserved and the observed log-likelihood is obtained by marginalization over all the label configurations:

$$ \mathcal{L}(\theta) = \log p(\mathbf{x}; \theta) = \log \left( \sum_{\mathbf{z} \in \mathcal{Z}, \mathbf{w} \in \mathcal{W}} p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \theta) \right). $$

Due to the double missing data structure $z$ for rows and $w$ for columns, neither the observed likelihood nor the E-step of the EM algorithm are tractable. Estimation can can nevertheless be performed either by numerical approximation, or by MCMC methods [see 8, 9].

## 2.2. Assumptions

We focus here on LBM where $\varphi$ belongs to a regular univariate exponential family set in canonical form:

$$ \varphi(x, \alpha) = b(x) \exp(\alpha x - \psi(\alpha)), $$

The canonical parameter $\alpha$ belongs to a space $\mathcal{A}$, so that $\varphi(\cdot, \alpha)$ is well defined for all $\alpha \in \mathcal{A}$. Classical properties of exponential families ensure that $\psi$ is convex, infinitely differentiable on $\hat{\mathcal{A}}$, and $(\psi')^{-1}$ is well defined on $\psi'(\hat{\mathcal{A}})$. When $X_\alpha \sim \varphi(\cdot, \alpha)$,

$$ \mathbb{E}[X_\alpha] = \psi'(\alpha) \text{ and } \mathbb{V}[X_\alpha] = \psi''(\alpha). $$
Notice that the definition of the exponential family used here relies on an exhaustive statistic that is X itself. This for a simple convenience. Family sets of the form $\varphi(x, \alpha) = b(x) \exp(\alpha t(x) - \psi(\alpha))$ can also be considered, all the further developments as Bernstein and concentration inequalities then concerning the exhaustive statistics $t(X)$.

Moreover, we make the following assumptions on the parameter space :

$H_1$: There exists a positive constant $c$, and a compact $C_\alpha$ such that

$$ \Theta \subset [c, 1-c]^g \times [c, 1-c]^m \times C_{\alpha}^{g \times m} \quad \text{with} \quad C_{\alpha} \subset \mathcal{A}. $$

$H_2$: The true parameter $\theta^* = (\pi^*, \rho^*, \alpha^*)$ lies in the relative interior of $\Theta$.

$H_3$: The mixture measure of LBM is identifiable: $\theta^*$ is identifiable up to a permutation of the row-labels and column-labels (see definition 2.7 of equivalent parameters).

The previous assumptions are standard. Notice that the following conditions are necessary for $H_3$ to hold:

$H_{3a}$: The map $\alpha \mapsto \varphi(\cdot, \alpha)$ is injective.

$H_{3b}$: Each row and each column of $\alpha^*$ is unique.

[9] gives sufficient conditions for the generic identifiability of the categorical LBM, i.e. except on a manifold set of null Lebesgue measure in $\Theta$ and this property is easily extended to the case of observations from a univariate exponential family. For binary SBM, [2] added the assumption $p \ge (\log n)/n$ on the parameter $p$ of the Bernoulli distribution to take into account sparsity.

Assumption $H_1$ ensures that the group proportions $\pi_k$ and $\rho_\ell$ are bounded away from 0 and 1 so that no group disappears when $n$ and $d$ go to infinity. It also ensures that $\alpha$ is bounded away from the boundaries of $\mathcal{A}$ and that there exists a positive value $\kappa > 0$, such that $[\alpha - \kappa, \alpha + \kappa] \subset \mathcal{A}$ for all parameters $\alpha$ of $C_\alpha$, which is essential to prove a uniform Bernstein inequality on the $X$.

Moreover, we define the quantity $\delta(\alpha)$ that captures the separation between row-groups or column-groups: low values of $\delta(\alpha)$ mean that two row-classes or two column-classes are very similar.

**Definition 2.3** (class distinctness). For $\theta = (\pi, \rho, \alpha) \in \Theta$. We define:

$$ \delta(\alpha) = \min \left\{ \min_{\ell \neq \ell'} \max_k \mathrm{KL}(\alpha_{k\ell}, \alpha_{k\ell'}), \min_{k \neq k'} \max_{\ell} \mathrm{KL}(\alpha_{k\ell}, \alpha_{k'\ell}) \right\} $$

with $\mathrm{KL}(\alpha, \alpha') = \mathbb{E}_\alpha[\log(\varphi(X, \alpha)/\varphi(X, \alpha'))] = \psi'(\alpha)(\alpha - \alpha') + \psi(\alpha') - \psi(\alpha)$ the Kullback divergence between $\varphi(\cdot, \alpha)$ and $\varphi(\cdot, \alpha')$.

*Remark 2.4.* Since $\alpha^*$ has distinct rows and distinct columns ($H_3$), $\delta(\alpha^*) > 0$.

*Remark 2.5.* These assumptions are satisfied for many distributions, including but not limited to:

* Bernoulli, when the proportion $p$ is bounded away from 0 and 1, or natural parameter $\alpha = \log(p/(1-p))$ bounded away from $\pm\infty$;
* Poisson, when the mean $\lambda$ is bounded away from 0 and $+\infty$, or natural parameter $\alpha = \log(\lambda)$ bounded away from $\pm\infty$;

* Gaussian with known variance when the mean $\mu$, which is also the natural parameter, is bounded away from $\pm\infty$.

In particular, the conditions stating that $\psi$ is twice differentiable and that $(\psi')^{-1}$ exists are equivalent to assuming that $X_\alpha$ has positive and finite variance for all values of $\alpha$ in the parameter space.

## 2.3. Model Symmetry

The LBM is a generalized mixture model and as such is subject to label switching. Moreover, the study of the asymptotics will involve the complete likelihood where symmetry properties on the parameter must be taken into account. We first recall the definition of a permutation in LBM, then define equivalence relationships for assignments and parameter, and discuss model symmetry.

**Definition 2.6 (permutation).** Let $s$ be a permutation on $\{1, \dots, g\}$ and $t$ a permutation on $\{1, \dots, m\}$. If $\mathbf{A}$ is a matrix with $g$ columns, we define $\mathbf{A}^s$ as the matrix obtained by permuting the columns of $\mathbf{A}$ according to $s$, i.e. for any row $i$ and column $k$ of $\mathbf{A}$, $A_{ik}^s = A_{is(k)}$. If $\mathbf{B}$ is a matrix with $m$ columns and $\mathbf{C}$ is a matrix with $g$ rows and $m$ columns, $\mathbf{B}^t$ and $\mathbf{C}^{s,t}$ are defined similarly:

$$ \mathbf{A}^s = (A_{is(k)})_{i,k} \quad \mathbf{B}^t = (B_{jt(l)})_{j,l} \quad \mathbf{C}^{s,t} = (C_{s(k)t(l)})_{k,l} $$

**Definition 2.7 (equivalence).** We define the following equivalence relationships:

* Two assignments $(\mathbf{z}, \mathbf{w})$ and $(\mathbf{z}', \mathbf{w}')$ are equivalent, denoted $\sim$, if they are equal up to label permutation, i.e. there exist two permutations $s$ and $t$ such that $\mathbf{z}' = \mathbf{z}^s$ and $\mathbf{w}' = \mathbf{w}^t$.

* Two parameters $\theta$ and $\theta'$ are equivalent, denoted $\sim$, if they are equal up to label permutation, i.e. there exist two permutations $s$ and $t$ such that $(\pi^s, \rho^t, \alpha^{s,t}) = (\pi', \rho', \alpha')$. This is label-switching.

* $(\theta, \mathbf{z}, \mathbf{w})$ and $(\theta', \mathbf{z}', \mathbf{w}')$ are equivalent, denoted $\sim$, if they are equal up to label permutation on $\alpha$, i.e. there exist two permutations, $s$ and $t$ such that $(\alpha^{s,t}, \mathbf{z}^s, \mathbf{w}^t) = (\alpha', \mathbf{z}', \mathbf{w}')$.

The last equivalence relationship is not concerned with $\pi$ and $\rho$. It is useful when dealing with the conditional likelihood $p(\mathbf{x}|\mathbf{z}, \mathbf{w}; \theta)$ which depends neither on $\pi$ nor $\rho$: in fact, if $(\theta, \mathbf{z}, \mathbf{w}) \sim (\theta', \mathbf{z}', \mathbf{w}')$, then for all $\mathbf{x}$, we have $p(\mathbf{x}|\mathbf{z}, \mathbf{w}; \theta) = p(\mathbf{x}|\mathbf{z}', \mathbf{w}'; \theta')$. Note also that $\mathbf{z} \sim \mathbf{z}^*$ (resp. $\mathbf{w} \sim \mathbf{w}^*$) if and only if there exists a permutation of the rows of the confusion matrix $\mathbb{R}_g(\mathbf{z})$ (resp. $\mathbb{R}_m(\mathbf{w})$) leading to a diagonal matrix.

**Definition 2.8 (symmetry).** We say that the parameter $\theta$ exhibits symmetry for the permutations $s, t$ if

$$ (\pi^s, \rho^t, \alpha^{s,t}) = (\pi, \rho, \alpha). $$
$\theta$ exhibits symmetry if it exhibits symmetry for any non trivial pair of permutations $(s, t)$. Finally the set of pairs $(s, t)$ for which $\theta$ exhibits symmetry is denoted $\text{Sym}(\theta)$.

**Remark 2.9.** The set of parameters that exhibit symmetry is a manifold of null Lebesgue measure in $\Theta$. This notion of symmetry is subtler than and different from label switching. To emphasize the difference between equivalence and symmetry, consider the following model: $\pi = (1/2, 1/2)$, $\rho = (1/3, 2/3)$ and $\alpha = \begin{pmatrix} \alpha_1 & \alpha_2 \\ \alpha_2 & \alpha_1 \end{pmatrix}$ with $\alpha_1 \neq \alpha_2$. The only permutations of interest here are $s = t = [1 2]$. Choose any $\mathbf{z}$ and $\mathbf{w}$. Because of label switching, we know that $p(\mathbf{x}, \mathbf{z}^s, \mathbf{w}^t; \theta^{s,t}) = p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \theta)$. $(\mathbf{z}^s, \mathbf{w}^t)$ and $(\mathbf{z}, \mathbf{w})$ have the same likelihood but under different parameters $\theta$ and $\theta^{s,t}$. If however, $\rho = (1/2, 1/2)$, then $(s, t) \in \text{Sym}(\theta)$ and $\theta^{s,t} = \theta$ so that $(\mathbf{z}, \mathbf{w})$ and $(\mathbf{z}^s, \mathbf{w}^t)$ have the same likelihood under the same parameter $\theta$. In particular, if $(\mathbf{z}, \mathbf{w})$ is a maximum-likelihood assignment under $\theta$, so is $(\mathbf{z}^s, \mathbf{w}^t)$. In other words, if $\theta$ exhibits symmetry, the maximum-likelihood assignment is not unique under the true model and there are at least $\#\text{Sym}(\theta)$ of them. This has important implications for the asymptotics of the observed likelihood ratio.

## 2.4. Distance and local assignments

We define the distance up to equivalence between two sets of assignments as follows:

**Definition 2.10 (distance).** The distance, up to equivalence, between configurations $\mathbf{z}$ and $\mathbf{z}^*$ is defined as

$$ \|\mathbf{z} - \mathbf{z}^*\|_{0,\sim} = \inf_{\mathbf{z}' \sim \mathbf{z}} \|\mathbf{z}' - \mathbf{z}^*\|_0 $$

where, for all matrix $\mathbf{z}$, $\|\cdot\|_0$ is the Hamming norm

$$ \|\mathbf{z}\|_0 = \sum_{i,k} 1\{z_{ik} \neq 0\}. $$

A similar definition is set for the distance between $\mathbf{w}$ and $\mathbf{w}^*$.

This allows us to define a neighborhood of radius $r$ in the assignment space, taking into account equivalent assignments classes.

**Definition 2.11 (Set of local assignments).** We denote $S(\mathbf{z}^*, \mathbf{w}^*, r)$ the set of configurations that have a representative (for $\sim$) within relative radius $r$ of $(\mathbf{z}^*, \mathbf{w}^*)$:

$$ S(\mathbf{z}^*, \mathbf{w}^*, r) = \{(\mathbf{z}, \mathbf{w}) : \|\mathbf{z} - \mathbf{z}^*\|_{0,\sim} \leq rn \text{ and } \|\mathbf{w} - \mathbf{w}^*\|_{0,\sim} \leq rd\} $$
### 3. Asymptotic properties in the complete data model

As stated in the introduction, we first study the asymptotic properties of the complete data model. Let $\hat{\theta}_c = (\hat{\pi}, \hat{\rho}, \hat{\alpha})$ be the MLE of $\theta$ in the complete data model, where the real assignments $\mathbf{z} = \mathbf{z}^*$ and $\mathbf{w} = \mathbf{w}^*$ are known. We can derive the following general estimates from Equation (2.1):

$$ \hat{\pi}_k(\mathbf{z}) = \frac{z_{+k}}{n} \qquad \hat{\rho}_\ell(\mathbf{w}) = \frac{w_{+\ell}}{d} \\ \hat{x}_{k\ell}(\mathbf{z}, \mathbf{w}) = \frac{\sum_{ij} x_{ij} z_{ik} w_{j\ell}}{z_{+k} w_{+\ell}} \qquad \hat{\alpha}_{k\ell} = \hat{\alpha}_{k\ell}(\mathbf{z}, \mathbf{w}) = (\psi')^{-1}(\hat{x}_{k\ell}(\mathbf{z}, \mathbf{w})) \qquad (3.1) $$

**Proposition 3.1.** The matrices $\Sigma_{\pi^*} = \text{Diag}(\pi^*) - \pi^* (\pi^*)^T$, $\Sigma_{\rho^*} = \text{Diag}(\rho^*) - \rho^* (\rho^*)^T$ are semi-definite positive, of rank $g-1$ and $m-1$, and $\hat{\pi}$ and $\hat{\rho}$ are asymptotically normal:

$$ \sqrt{n}(\hat{\pi}(\mathbf{z}^*) - \pi^*) \xrightarrow[n \to \infty]{D} \mathcal{N}(0, \Sigma_{\pi^*}) \quad \text{and} \quad \sqrt{d}(\hat{\rho}(\mathbf{w}^*) - \rho^*) \xrightarrow[d \to \infty]{D} \mathcal{N}(0, \Sigma_{\rho^*}) \qquad (3.2) $$

Similarly, let $V(\alpha^*)$ be the matrix defined by $[V(\alpha^*)]_{k\ell} = 1/\psi''(\alpha_{k\ell}^*)$ and $\Sigma_{\alpha^*} = \text{Diag}^{-1}(\pi^*)V(\alpha^*)\text{Diag}^{-1}(\rho^*)$. Then:

$$ \sqrt{nd}(\hat{\alpha}_{k\ell}(\mathbf{z}^*, \mathbf{w}^*) - \alpha_{k\ell}^*) \xrightarrow[n,d \to \infty]{D} \mathcal{N}(0, \Sigma_{\alpha^*, k\ell}) \quad \text{for all } k, \ell $$

and the components $\hat{\alpha}_{k,\ell}$ are independent.

*Proof:* Since $\hat{\pi}(\mathbf{z}^*) = (\hat{\pi}_1(\mathbf{z}^*), \dots, \hat{\pi}_g(\mathbf{z}^*))$ (resp. $\hat{\rho}(\mathbf{w}^*)$) is the sample mean of $n$ (resp. $d$) i.i.d. multinomial random variables with parameters $1$ and $\pi^*$ (resp. $\rho^*$), a simple application of the central limit theorem (CLT) gives:

$$ \Sigma_{\pi^*,kk'} = \begin{cases} \pi_k^*(1 - \pi_k^*) & \text{if } k = k' \\ -\pi_k^*\pi_{k'}^* & \text{if } k \neq k' \end{cases} \quad \text{and} \quad \Sigma_{\rho^*,\ell\ell'} = \begin{cases} \rho_\ell^*(1 - \rho_\ell^*) & \text{if } \ell = \ell' \\ -\rho_\ell^*\rho_{\ell'}^* & \text{if } \ell \neq \ell' \end{cases} $$

which proves Equation (3.2) where $\Sigma_{\pi^*}$ and $\Sigma_{\rho^*}$ are semi-definite positive of rank $g-1$ and $m-1$.

Similarly, $\psi'(\hat{\alpha}_{k\ell}(\mathbf{z}^*, \mathbf{w}^*))$ is the average of $z_{+k}^* w_{+\ell}^* = nd\hat{\pi}_k(\mathbf{z}^*) \hat{\rho}_\ell(\mathbf{w}^*)$ i.i.d. random variables with mean $\psi'(\alpha_{k\ell}^*)$ and variance $\psi''(\alpha_{k\ell}^*)$. $nd\hat{\pi}_k(\mathbf{z}^*) \hat{\rho}_\ell(\mathbf{w}^*)$ is itself random but $\hat{\pi}_k(\mathbf{z}^*) \hat{\rho}_\ell(\mathbf{w}^*) \xrightarrow[n,d \to +\infty]{D} \pi_k^* \rho_\ell^*$ almost surely. Therefore, by Slutsky's lemma and the CLT for random sums of random variables [12], we have:

$$ \sqrt{nd\hat{\pi}_k^* \rho_\ell^*} (\psi'(\hat{\alpha}_{k\ell}(\mathbf{z}^*, \mathbf{w}^*)) - \psi'(\alpha_{k\ell}^*)) = \sqrt{nd\hat{\pi}_k^* \rho_\ell^*} \left( \frac{\sum_{ij} X_{ij} z_{ik}^* w_{j\ell}^*}{nd\hat{\pi}_k(\mathbf{z}^*) \hat{\rho}_\ell(\mathbf{w}^*)} - \psi'(\alpha_{k\ell}^*) \right) \xrightarrow[n,d \to +\infty]{D} \mathcal{N}(0, \psi''(\alpha_{k\ell}^*)) $$
The differentiability of $(\psi')^{-1}$ and the delta method then gives:

$$ \sqrt{nd}(\hat{\alpha}_{k\ell}(\mathbf{z}^*, \mathbf{w}^*) - \alpha_{k\ell}^*) \xrightarrow[n,d \to +\infty]{D} \mathcal{N}\left(0, \frac{1}{\pi_k^* \rho_\ell^* \psi''(\alpha_{k\ell}^*)}\right) $$

and the independence results from the independence of $\hat{\alpha}_{k\ell}(\mathbf{z}^*, \mathbf{w}^*)$ and $\hat{\alpha}_{k'\ell'}(\mathbf{z}^*, \mathbf{w}^*)$ as soon as $k \neq k'$ or $\ell \neq \ell'$, as they involve different sets of independent variables.

□

Moreover, the complete model is locally asymptotically normal (LAN), as stated in the following proposition. Note that the unusual condition for $s$ and $t$ arises from the constraints $\pi^T \mathbf{1}_g = \rho^T \mathbf{1}_m = 1$, where $\mathbf{1}_g$ is the vector of size $g$ filled with 1, which must be satisfied even after perturbing $\pi$ (resp. $\rho$) with $s$ (resp. $t$).

**Proposition 3.2 (Local asymptotic normality).** Let $\mathcal{L}_c^*$ the map defined by $\theta = (\pi, \rho, \alpha) \mapsto \log p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta)$ and note $I_{\pi^*} = \text{Diag}^{-1}(\pi^*)$, $I_{\rho^*} = \text{Diag}^{-1}(\rho^*)$ and $I_{\alpha^*}$ the component-wise inverse of $\Sigma_{\alpha^*}$. For any $s$, $t$ and $u$ in a compact set, such that $t^T \mathbf{1}_g = 0$ and $s^T \mathbf{1}_m = 0$, we have:

$$ \mathcal{L}_c^* \left( \pi^* + \frac{s}{\sqrt{n}}, \rho^* + \frac{t}{\sqrt{d}}, \alpha^* + \frac{u}{\sqrt{nd}} \right) = \mathcal{L}_c^* (\theta^*) + s^T \mathbf{Y}_{\pi^*} + t^T \mathbf{Y}_{\rho^*} + \text{Tr}(u^T \mathbf{Y}_{\alpha^*}) \\ - \left( \frac{1}{2} s^T I_{\pi^*} s + \frac{1}{2} t^T I_{\rho^*} t + \frac{1}{2} \text{Tr}[(u \odot u)^T I_{\alpha^*}] \right) \\ + o_P(1) $$

where $\odot$ denotes the Hadamard product of two matrices (element-wise product), $\mathbf{Y}_{\pi^*}$, $\mathbf{Y}_{\rho^*}$ are asymptotically centered Gaussian vectors of sizes $g$ and $m$ with respective variance matrices $I(\pi^*)$ and $I(\rho^*)$ and $\mathbf{Y}_{\alpha^*}$ is a random matrix of size $g \times m$ with independent Gaussian components $Y_{\alpha^*,kl} \sim \mathcal{N}(0, I_{\alpha^*,kl})$.

**Proof.**

By Taylor expansion, and with the condition $s^T \mathbf{1}_g = t^T \mathbf{1}_m = 0$

$$ \mathcal{L}_c^* \left( \pi^* + \frac{s}{\sqrt{n}}, \rho^* + \frac{t}{\sqrt{d}}, \alpha^* + \frac{u}{\sqrt{nd}} \right) = \mathcal{L}_c^* (\theta^*) + \frac{1}{\sqrt{n}} s^T \nabla \mathcal{L}_{c\pi}^* (\theta^*) + \frac{1}{\sqrt{d}} t^T \nabla \mathcal{L}_{c\rho}^* (\theta^*) + \frac{1}{\sqrt{nd}} \text{Tr} (u^T \nabla \mathcal{L}_{c\alpha}^* (\theta^*)) \\ + \frac{1}{n} s^T \mathbf{H}_\pi (\theta^*) s + \frac{1}{d} t^T \mathbf{H}_\rho (\theta^*) t + \frac{1}{nd} \text{Tr}((u \odot u)^T \mathbf{H}_\alpha (\theta^*)) + o_P(1) $$

where $\nabla \mathcal{L}_{c\pi}^* (\theta^*)$, $\nabla \mathcal{L}_{c\rho}^* (\theta^*)$ and $\nabla \mathcal{L}_{c\alpha}^* (\theta^*)$ denote the respective components of the gradient of $\mathcal{L}_c^*$ evaluated at $\theta^*$ and $\mathbf{H}_\pi$, $\mathbf{H}_\rho$ and $\mathbf{H}_\alpha$ denotes the conditional hessian of $\mathcal{L}_c^*$ evaluated at $\theta^*$. By inspection, $\mathbf{H}_\pi/n$, $\mathbf{H}_\rho/d$ and $\mathbf{H}_\alpha/nd$
converge in probability to constant matrices and the random vectors $\nabla \mathcal{L}_{c\pi}^*(\theta^*) / \sqrt{n}$, $\nabla \mathcal{L}_{c\rho}^*(\theta^*) / \sqrt{d}$ and $\nabla \mathcal{L}_{c\alpha}^*(\theta^*) / \sqrt{nd}$ converge in distribution to Gaussian vectors by the central limit theorem.

<figure><img src="image_4.png" /></figure>

# 4. Profile Likelihood

Our main result compares the observed likelihood ratio $p(\mathbf{x}; \theta)/p(\mathbf{x}; \theta^*)$ with the complete likelihood $p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta)/p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*)$. To study the behavior of these likelihoods, we shall work conditionally to the true configurations $(\mathbf{z}^*, \mathbf{w}^*)$ that have enough observations in each row or column group. We therefore define in section 4.1 so called regular configurations and prove that they occur with high probability. We then introduce in section 4.2 conditional and profile likelihood ratios and state some of their properties.

## 4.1. Regular assignments

**Definition 4.1** (c-regular assignments). Let $\mathbf{z} \in \mathcal{Z}$ and $\mathbf{w} \in \mathcal{W}$. For any $c > 0$, we say that $\mathbf{z}$ and $\mathbf{w}$ are c-regular if

$$ \min_k z_{+k} \ge cn \quad \text{and} \quad \min_l w_{+l} \ge cd. $$

In regular configurations, each row-group for example has $\Omega(n)$ members, where $u_n = \Omega(n)$ if there exists two constant $a, b > 0$ such that for $n$ enough large $an \le u_n \le bn$. c/2-regular assignments, with $c$ defined in Assumption $H_1$, have high $\mathbb{P}_{\theta^*}$-probability in the space of all assignments, uniformly over all $\theta^* \in \Theta$, as stated in Proposition 4.2.

**Proposition 4.2.** Define $\mathcal{Z}_1$ and $\mathcal{W}_1$ as the subsets of $\mathcal{Z}$ and $\mathcal{W}$ made of c/2-regular assignments, with $c$ defined in assumption $H_1$. Denote $\Omega_1$ the event $\{(\mathbf{z}^*, \mathbf{w}^*) \in \mathcal{Z}_1 \times \mathcal{W}_1\}$, then:

$$ \mathbb{P}_{\theta^*}(\bar{\Omega}_1) \le g \exp\left(-\frac{nc^2}{2}\right) + m \exp\left(-\frac{dc^2}{2}\right). $$

Each $z_{+k}$ is a sum of $n$ i.i.d Bernoulli random variables with parameter $\pi_k \ge \pi_{\min} \ge c$. The proof is straightforward and stems from a simple Hoeffding bound

$$ \mathbb{P}_{\theta^*} \left( z_{+k} \le n \frac{c}{2} \right) \le \mathbb{P}_{\theta^*} \left( z_{+k} \le n \frac{\pi_k}{2} \right) \le \exp \left( -2n \left( \frac{\pi_k}{2} \right)^2 \right) \le \exp \left( -\frac{nc^2}{2} \right) $$

and a union bound over $g$ values of $k$, with similar approach for $w_{+l}$.
## 4.2. Conditional and profile log-likelihoods

Introducing the conditional log-likelihood ratio

$$F_{nd}(\theta, \mathbf{z}, \mathbf{w}) = \log \frac{p(\mathbf{x}|\mathbf{z}, \mathbf{w}; \theta)}{p(\mathbf{x}|\mathbf{z}^*, \mathbf{w}^*; \theta^*)},$$

the complete likelihood can be written as follows

$$p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \theta) = p(\mathbf{z}, \mathbf{w}; \theta)p(\mathbf{x}|\mathbf{z}^*, \mathbf{w}^*; \theta^*) \exp(F_{nd}(\theta, \mathbf{z}, \mathbf{w})).$$

The study of $F_{nd}$ will be of crucial importance, as well as its maximum over $\Theta$. After some definitions, we examine some useful properties.

**Definition 4.3.** The conditional expectation $G$ of $F_{nd}$ is defined as:

$$G(\theta, \mathbf{z}, \mathbf{w}) = \mathbb{E}_{\theta} \left[ \log \frac{p(\mathbf{x}|\mathbf{z}, \mathbf{w}; \theta)}{p(\mathbf{x}|\mathbf{z}^*, \mathbf{w}^*; \theta^*)} \middle| \mathbf{z}^*, \mathbf{w}^* \right] = \mathbb{E}_{\theta} \left[ F_{nd}(\theta, \mathbf{z}, \mathbf{w}) \middle| \mathbf{z}^*, \mathbf{w}^* \right]$$

Moreover, the profile log-likelihood ratio $\Lambda$ and its expectation $\tilde{\Lambda}$ are defined as:

$$\Lambda(\mathbf{z}, \mathbf{w}) = \max_{\theta} F_{nd}(\theta, \mathbf{z}, \mathbf{w})$$

$$\tilde{\Lambda}(\mathbf{z}, \mathbf{w}) = \max_{\theta} G(\theta, \mathbf{z}, \mathbf{w}).$$

**Remark 4.4.** As $F_{nd}$ and $G$ only depend on $\theta$ through $\alpha$, we will sometimes replace $\theta$ with $\alpha$ in the expressions of $F_{nd}$ and $G$. Replacing $F_{n,d}$ and $G$ by their profiled version $\Lambda$ and $\tilde{\Lambda}$ allows us to get rid of the continuous argument $\alpha$ of $F_{nd}$ and to rely instead only on discrete contrasts $\Lambda$ and $\tilde{\Lambda}$.

Now, Proposition 4.5 characterizes which values of $\alpha$ maximize $F_{nd}$ and $G$ to reach $\Lambda$ and $\tilde{\Lambda}$. Propositions 4.6 and 4.7 in turn describes properties of $G$ and $\tilde{\Lambda}$ relative to $(\mathbf{z}, \mathbf{w})$.

**Proposition 4.5** (maximum of $G$ and $\tilde{\Lambda}$ in $\theta$). Let $\hat{\theta}_c = (\hat{x}_{k\ell}(\mathbf{z}, \mathbf{w}), \hat{\pi}_k(\mathbf{z}), \hat{\rho}_\ell(\mathbf{w}))$ be the maximum likelihood estimator of the complete model, as defined in Equation 3.1. Conditionally on $\mathbf{z}^*, \mathbf{w}^*$, define the following quantities:

$$S^* = (S_{k\ell}^*)_{k\ell} = (\psi'(\alpha_{k\ell}^*))_{k\ell}$$

$$\bar{x}_{k\ell}(\mathbf{z}, \mathbf{w}) = \mathbb{E}_{\theta^*} [\hat{x}_{k\ell}(\mathbf{z}, \mathbf{w}) | \mathbf{z}^*, \mathbf{w}^*] = \frac{[\mathbb{R}_g(\mathbf{z})^T S^* \mathbb{R}_m(\mathbf{w})]_{k\ell}}{\hat{\pi}_k(\mathbf{z}) \hat{\rho}_\ell(\mathbf{w})} \quad (4.1)$$

with $\bar{x}_{k\ell}(\mathbf{z}, \mathbf{w}) = 0$ for $\mathbf{z}$ and $\mathbf{w}$ such that $\hat{\pi}_k(\mathbf{z}) = 0$ or $\hat{\rho}_\ell(\mathbf{w}) = 0$.

Then $F_{nd}(\theta, \mathbf{z}, \mathbf{w})$ (resp. $G(\theta, \mathbf{z}, \mathbf{w})$) is maximum in $\alpha$ for $\alpha = \hat{\alpha}(\mathbf{z}, \mathbf{w})$ (resp. $\bar{\alpha}(\mathbf{z}, \mathbf{w})$) defined by:

$$\hat{\alpha}(\mathbf{z}, \mathbf{w})_{k\ell} = (\psi')^{-1}(\hat{x}_{k\ell}(\mathbf{z}, \mathbf{w})) \quad \text{and} \quad \bar{\alpha}(\mathbf{z}, \mathbf{w})_{k\ell} = (\psi')^{-1}(\bar{x}_{k\ell}(\mathbf{z}, \mathbf{w})).$$

Hence,

$$\Lambda(\mathbf{z}, \mathbf{w}) = F_{nd}(\hat{\alpha}(\mathbf{z}, \mathbf{w}), \mathbf{z}, \mathbf{w})$$

$$\tilde{\Lambda}(\mathbf{z}, \mathbf{w}) = G(\bar{\alpha}(\mathbf{z}, \mathbf{w}), \mathbf{z}, \mathbf{w})$$
Note that although $\bar{x}_{k\ell} = \mathbb{E}_{\theta^*}[\hat{x}_{k\ell}|\mathbf{z}^*, \mathbf{w}^*]$, in general $\bar{\alpha}_{k\ell} \neq \mathbb{E}_{\theta^*}[\hat{\alpha}_{k\ell}|\mathbf{z}^*, \mathbf{w}^*]$ by non linearity of $(\psi')^{-1}$. Nevertheless, since $(\psi')^{-1}$ is Lipschitz over compact subsets of $\psi'(\tilde{\Lambda})$, with high probability, $|\bar{\alpha}_{k\ell} - \hat{\alpha}_{k\ell}|$ and $|\hat{x}_{k\ell} - \bar{x}_{k\ell}|$ are of the same order of magnitude.

**Proposition 4.6** (maximum of $G$ and $\tilde{\Lambda}$ in $(\theta, \mathbf{z}, \mathbf{w})$). Let KL$(\alpha, \alpha') = \psi'(\alpha)(\alpha - \alpha') + \psi(\alpha') - \psi(\alpha)$ be the Kullback divergence between $\varphi(., \alpha)$ and $\varphi(., \alpha')$ then:

$$G(\theta, \mathbf{z}, \mathbf{w}) = -nd \sum_{k,k'} \sum_{\ell,\ell'} \mathbb{R}_g(\mathbf{z})_{k,k'} \mathbb{R}_m(\mathbf{w})_{\ell,\ell'} \mathrm{KL}(\alpha_{k\ell}^*, \alpha_{k'\ell'}) \le 0. \quad (4.2)$$

Conditionally on the set $\Omega_1$ of regular assignments and for $n, d > 2/c$,

(i) $G$ is maximized at $(\alpha^*, \mathbf{z}^*, \mathbf{w}^*)$ and its equivalence class.

(ii) $\tilde{\Lambda}$ is maximized at $(\mathbf{z}^*, \mathbf{w}^*)$ and its equivalence class and $\tilde{\Lambda}(\mathbf{z}^*, \mathbf{w}^*) = 0$.

Moreover, the maximum of $\tilde{\Lambda}$ in $(\mathbf{z}^*, \mathbf{w}^*)$ is well separated, in the sense that there exists a positive gap between $\tilde{\Lambda}(\mathbf{z}^*, \mathbf{w}^*)$ and any other $\tilde{\Lambda}(\mathbf{z}, \mathbf{w})$ for $(\mathbf{z}, \mathbf{w})$ in a close neighborhood of $(\mathbf{z}^*, \mathbf{w}^*)$, as stated in the following proposition:

**Proposition 4.7** (Separability for $\tilde{\Lambda}$). Conditionally upon $\Omega_1$, there exists a positive constant $C$ such that for all $(\mathbf{z}, \mathbf{w}) \in S(\mathbf{z}^*, \mathbf{w}^*, C)$:

$$\tilde{\Lambda}(\mathbf{z}, \mathbf{w}) \le -\frac{c\delta(\alpha^*)}{4} (d\|\mathbf{z} - \mathbf{z}^*\|_{0,\sim} + n\|\mathbf{w} - \mathbf{w}^*\|_{0,\sim}) \quad (4.3)$$

Moreover, there exists a positive constant $B(C)$ such that for all $(\mathbf{z}, \mathbf{w}) \notin S(\mathbf{z}^*, \mathbf{w}^*, C)$

$$\tilde{\Lambda}(\mathbf{z}, \mathbf{w}) \le -B(C) nd \quad (4.4)$$

The proofs of these propositions are reported in Appendix A. Proof of Proposition 4.5 follows from a straightforward calculation, proof of Proposition 4.6 uses the technical Lemma D.1 to characterize the maximum of $G$ and proof of Proposition 4.7 uses regularity properties of the gradient of $\tilde{\Lambda}$ to control its behavior near its maximum.

## 5. Main Result

Our main result matches the asymptotics of complete and observed likelihoods and is the key to prove the consistency of maximum likelihood and variational estimators. It is set under the assumptions described in section 2.2 and the following asymptotics for the number of rows $n$ and columns $d$:

$$(H_4): \log(d)/n \to 0 \text{ and } \log(n)/d \to 0.$$

**Theorem 5.1** (complete-observed). Let $\mathbf{x}$ be a matrix of $n \times d$ observations of a LBM with true parameter $\theta^* = (\pi^*, \rho^*, \alpha^*)$ where the number of row-groups $g$ and column-groups $m$ are known, which conditional distribution belongs to a regular univariate exponential family. The true random and unobserved
assignations for rows and columns are denoted $\mathbf{z}^*$ and $\mathbf{w}^*$ respectively. Define $\#\text{Sym}(\theta)$ as the number of pairs of permutations $(s,t)$ for which $\theta$ exhibits symmetry.

If assumptions $H_1$ to $H_4$ are fulfilled, then, the observed likelihood ratio behaves like the complete likelihood ratio, up to a bounded multiplicative factor:

$$ \frac{p(\mathbf{x}; \boldsymbol{\theta})}{p(\mathbf{x}; \boldsymbol{\theta}^*)} = \frac{\#\text{Sym}(\boldsymbol{\theta})}{\#\text{Sym}(\boldsymbol{\theta}^*)} \max_{\boldsymbol{\theta}' \sim \boldsymbol{\theta}} \frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}')}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*)} (1 + o_P(1)) + o_P(1) $$

where both $o_P$ are uniform over all $\boldsymbol{\theta} \in \Theta$.

The maximum over all $\boldsymbol{\theta}'$ that are equivalent to $\boldsymbol{\theta}$ stems from the fact that because of label-switching, $\boldsymbol{\theta}$ is only identifiable up to its $\sim$-equivalence class from the observed likelihood, whereas it is completely identifiable from the complete likelihood as in this latter case, the labels are known. The terms $\#\text{Sym}$ are needed to take into account cases where $\boldsymbol{\theta}$ exhibits symmetry. These were omitted by [2] for SBM, although they are also needed in this case, see remark 5.3. When no $\boldsymbol{\theta} \in \Theta$ exhibits symmetry, the following corollary is immediately deduced:

**Corollary 5.2.** If $\Theta$ contains only parameters that do not exhibit symmetry:

$$ \frac{p(\mathbf{x}; \boldsymbol{\theta})}{p(\mathbf{x}; \boldsymbol{\theta}^*)} = \max_{\boldsymbol{\theta}' \sim \boldsymbol{\theta}} \frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}')}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*)} (1 + o_P(1)) + o_P(1) $$

where the $o_P$ is uniform over all $\boldsymbol{\theta} \in \Theta$.

**General sketch of the proof.** The proof relies on the following decomposition of the observed likelihood:

$$ p(\mathbf{x}; \boldsymbol{\theta}) = \sum_{(\mathbf{z}, \mathbf{w})} p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \boldsymbol{\theta}) = \sum_{(\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)} p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \boldsymbol{\theta}) + \sum_{(\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)} p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \boldsymbol{\theta}). $$

where the second term shall be proved to be asymptotically negligible. Its control stems from the study of the conditional log-likelihood $F_{nd}$, see Equation 4.2. In fact, the contribution of configurations that are not equivalent to $(\mathbf{z}^*, \mathbf{w}^*)$ leads itself to the study of a global control, and a sharper local control of $F_{nd}$. Hence, the proof relies on the examination of the asymptotic behavior of $F_{nd}$ on three types of configurations that partition $\mathcal{Z} \times \mathcal{W}$:

1. **global control** for assignations $(\mathbf{z}, \mathbf{w})$ sufficiently far from $(\mathbf{z}^*, \mathbf{w}^*)$, i.e. such that $\tilde{\Lambda}(\mathbf{z}, \mathbf{w})$ is of order $\Omega(-nd)$. Proposition 5.5 gives a large deviation result for $F_{nd} - \tilde{\Lambda}(\mathbf{z}, \mathbf{w})$ to prove that $F_{nd}$ is also of order $-\Omega_P(nd)$. A key point will be the use of Proposition C.4, establishing a specific concentration inequality for sub-exponential variables. In turn, those assignments contribute as a $o_P(p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*))$ to the sum (Proposition 5.6).

2. **local control**: a small deviation result (Proposition 5.7) is needed to show that the combined contribution of assignments close to but not equivalent to $(\mathbf{z}^*, \mathbf{w}^*)$ is also a $o_P(p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*))$ (Proposition 5.8).
3. equivalent assignments: Proposition 5.9 examines which of the remaining assignments, all equivalent to $(\mathbf{z}^*, \mathbf{w}^*)$, contribute to the sum.

Once these propositions proved, the proof is straightforward, as can be seen below. They are in turn carefully presented and discussed in dedicated sub-sections as they represents the core arguments and their proofs are themselves postponed to Appendix B for more readability.

**Proof.**

We work conditionally to $\Omega_1$, defined in Proposition 4.2, i.e., the high probability event that $(\mathbf{z}^*, \mathbf{w}^*)$ is a c/2-regular assignment. We choose $(\mathbf{z}^*, \mathbf{w}^*) \in \mathcal{Z}_1 \times W_1$ and a sequence $t_{nd}$ decreasing to 0 but satisfying $t_{nd}^2 \gg \frac{n+d}{nd}$. This is possible when $n \to \infty$ and $d \to \infty$, and for example with Assumption (H$_4$). We write:

$$p(\mathbf{x}; \boldsymbol{\theta}) = \sum_{(\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)} p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \boldsymbol{\theta}) + \sum_{(\mathbf{z}, \mathbf{w}) \notin S(\mathbf{z}^*, \mathbf{w}^*, t_{nd})} p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \boldsymbol{\theta}) + \sum_{\substack{(\mathbf{z}, \mathbf{w}) \in S(\mathbf{z}^*, \mathbf{w}^*, t_{nd}) \\ (\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)}} p(\mathbf{z}, \mathbf{w}, \mathbf{x}; \boldsymbol{\theta})$$

According to Proposition 5.6, conditionally to $\Omega_1$ and for $n, d$ large enough that $2\sqrt{2nd}t_{nd} \ge gm$, the contribution of far away assignments is

$$\sup_{\boldsymbol{\theta} \in \Theta} \sum_{(\mathbf{z}, \mathbf{w}) \notin S(\mathbf{z}^*, \mathbf{w}^*, t_{nd})} p(\mathbf{z}, \mathbf{w}, \mathbf{x}; \boldsymbol{\theta}) = o_P(p(\mathbf{z}^*, \mathbf{w}^*, \mathbf{x}; \boldsymbol{\theta}^*)).$$

Using the separability of $\tilde{\Lambda}$ and Assumption (H$_4$), Proposition 5.8 ensures the existence of $C > 0$ such that:

$$\sup_{\boldsymbol{\theta} \in \Theta} \sum_{\substack{(\mathbf{z}, \mathbf{w}) \in S(\mathbf{z}^*, \mathbf{w}^*, C) \\ (\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)}} p(\mathbf{z}, \mathbf{w}, \mathbf{x}; \boldsymbol{\theta}) = o_P(p(\mathbf{z}^*, \mathbf{w}^*, \mathbf{x}; \boldsymbol{\theta}^*))$$

Since $t_{nd}$ decreases to 0, Proposition 5.8 can be applied for the local configurations belonging to $S(t_{nd})$, for $n, d$ large enough. Therefore the observed likelihood ratio reduces to:

$$\frac{p(\mathbf{x}; \boldsymbol{\theta})}{p(\mathbf{x}; \boldsymbol{\theta}^*)} = \frac{\sum_{(\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)} p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \boldsymbol{\theta}) + p(\mathbf{x}; \mathbf{z}^*, \mathbf{w}^*, \boldsymbol{\theta}^*) o_P(1)}{\sum_{(\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)} p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \boldsymbol{\theta}^*) + p(\mathbf{x}; \mathbf{z}^*, \mathbf{w}^*, \boldsymbol{\theta}^*) o_P(1)}$$

Proposition 5.9 deals with equivalence and symmetry and allows us to conclude

$$\frac{p(\mathbf{x}; \boldsymbol{\theta})}{p(\mathbf{x}; \boldsymbol{\theta}^*)} = \frac{\# \mathrm{Sym}(\boldsymbol{\theta})}{\# \mathrm{Sym}(\boldsymbol{\theta}^*)} \max_{\boldsymbol{\theta}' \sim \boldsymbol{\theta}} \frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}')}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*)} (1 + o_P(1)) + o_P(1).$$
<figure><img src="image_5.png" /></figure>

Remark 5.3. As already pointed out, if $\theta$ exhibits symmetry, the maximum likelihood assignment is not unique under $\theta$, and $\# \text{Sym}(\theta)$ terms contribute with the same weight. This was not taken into account by [2], and it is interesting to see why it should be also present for SBM. Recall that SBM has only one set of labels $z$. The proof relies on the the decomposition

$$p(\mathbf{x}; \boldsymbol{\theta}) = \sum_{\mathbf{z}} p(\mathbf{x}, \mathbf{z}; \boldsymbol{\theta}) = \sum_{\mathbf{z}' \sim \mathbf{z}^*} p(\mathbf{x}, \mathbf{z}'; \boldsymbol{\theta}) + \sum_{\mathbf{z}' \sim \mathbf{z}^*} p(\mathbf{x}, \mathbf{z}'; \boldsymbol{\theta})$$

where the second term of the sum is neglectible compared to the first term. Now, $\mathbf{z}' \sim \mathbf{z}^*$ means that there exists a permutation $t: [g] \to [g]$ such that $\mathbf{z}' = \mathbf{z}^t$ and $p(\mathbf{x}, \mathbf{z}^t; \boldsymbol{\theta}) = p(\mathbf{x}, \mathbf{z}; \boldsymbol{\theta}^t)$. The first term is written on Page 1941, Equation (25) in [2] as

$$\sum_{\mathbf{z}' \sim \mathbf{z}^*} p(\mathbf{x}, \mathbf{z}'; \boldsymbol{\theta}) = \sum_{\boldsymbol{\theta}' \sim \boldsymbol{\theta}} p(\mathbf{x}, \mathbf{z}^*; \boldsymbol{\theta}') = (1 + o(1)) \max_{\boldsymbol{\theta}' \sim \boldsymbol{\theta}} p(\mathbf{x}, \mathbf{z}^*; \boldsymbol{\theta}')$$

However, the first equality is not always correct. Actually, we have

$$\sum_{\mathbf{z}' \sim \mathbf{z}^*} p(\mathbf{x}, \mathbf{z}'; \boldsymbol{\theta}) = \sum_{t: [g] \to [g]} p(\mathbf{x}, \mathbf{z}^{*,t}; \boldsymbol{\theta}) = \sum_{t: [g] \to [g]} p(\mathbf{x}, \mathbf{z}^*; \boldsymbol{\theta}^t)$$

Take a special case of symmetry where $\pi = (1/g, \dots, 1/g)$ and $\alpha = (p-q)I_g + q1_g 1_g^T$. Then we have $\boldsymbol{\theta}^t = \boldsymbol{\theta}$ for all $t$. Thus,

$$\sum_{\mathbf{z}' \sim \mathbf{z}^*} p(\mathbf{x}, \mathbf{z}'; \boldsymbol{\theta}) = g! p(\mathbf{x}, \mathbf{z}^*; \boldsymbol{\theta}).$$

Even for the SBM, we thus have generally:

$$\sum_{\mathbf{z}' \sim \mathbf{z}^*} p(\mathbf{x}, \mathbf{z}'; \boldsymbol{\theta}) = (1 + o(1)) \# \text{Sym}(\boldsymbol{\theta}) \max_{\boldsymbol{\theta}' \sim \boldsymbol{\theta}} p(\mathbf{x}, \mathbf{z}^*; \boldsymbol{\theta}')$$

## 5.1. Global Control

A large deviation inequality for configurations $(\mathbf{z}, \mathbf{w})$ far from $(\mathbf{z}^*, \mathbf{w}^*)$ is built and used to prove that far away configurations make a small contribution to $p(\mathbf{x}; \boldsymbol{\theta})$. Since we restricted $\alpha$ in a bounded subset of $\hat{\mathcal{A}}$, there exists two positive values $M_\alpha$ and $\kappa$ such that $C_\alpha + (-\kappa, \kappa) \subset [-M_\alpha, M_\alpha] \subset \hat{\mathcal{A}}$. Moreover, the variance of $X_\alpha$ is bounded away from 0 and $+\infty$:

$$\sup_{\alpha \in [-M_\alpha, M_\alpha]} \mathbb{V}(X_\alpha) = \bar{\sigma}^2 < +\infty \quad \text{and} \quad \inf_{\alpha \in [-M_\alpha, M_\alpha]} \mathbb{V}(X_\alpha) = \underline{\sigma}^2 > 0.$$

**Proposition 5.4.** With the previous notations, if $\alpha \in C_\alpha$ and $X_\alpha \sim \varphi(., \alpha)$, then $X_\alpha$ is sub-exponential with parameters $(\bar{\sigma}^2, \kappa^{-1})$.
The latter proposition is a direct consequence of the definition of sub-exponential variables, see Appendix C.

**Proposition 5.5** (large deviations of $F_{nd}$). Let $\text{Diam}(\Theta) = \sup_{\theta, \theta'} \|\theta - \theta'\|_{\infty}$. For all $\varepsilon_{n,d} \le \kappa\bar{\sigma}$ and $n, d$

$$ \begin{aligned} & \Delta_{nd}^{1}(\varepsilon_{nd}) \\ & = \mathbb{P} \left( \sup_{\theta, \mathbf{z}, \mathbf{w}} \left\{ F_{nd}(\theta, \mathbf{z}, \mathbf{w}) - \tilde{\Lambda}(\mathbf{z}, \mathbf{w}) \right\} \ge \bar{\sigma}nd \text{Diam}(\Theta) \sqrt{2\varepsilon_{nd}} \left[ 1 + \frac{gm}{2\sqrt{2nd\varepsilon_{nd}}} \right] \right) \\ & \le g^{n}m^{d} \exp \left( -\frac{nd\varepsilon_{nd}^{2}}{2} \right) \quad (5.1) \end{aligned} $$

In particular, if $n$ and $d$ are large enough that $2\sqrt{2nd\varepsilon_{nd}} \ge gm$, the previous inequality ensures that with high probability, $F_{nd}(\theta, \mathbf{z}, \mathbf{w}) - \tilde{\Lambda}(\mathbf{z}, \mathbf{w})$ is no greater than $\bar{\sigma}nd \text{Diam}(\Theta)4\sqrt{2\varepsilon_{nd}}$.

The concentration inequality used in [2] to prove an analog result for SBM is not sufficient here, as it can be used only for upper-bounded observations, which is obviously not the case for all exponential families. We instead develop a Bernstein-type inequality for sub-exponential variables (Proposition C.4) to upper bound $F_{nd}(\theta, \mathbf{z}, \mathbf{w}) - \tilde{\Lambda}(\mathbf{z}, \mathbf{w})$. Proposition 5.5 relies heavily on this Bernstein inequality. A straightforward consequence of this deviation bound is that the combined contribution of assignments far away from $(\mathbf{z}^*, \mathbf{w}^*)$ to the sum is negligible, assuming that the numbers $n$ of rows and $d$ of columns grow at commensurate rates, as stated in the following proposition:

**Proposition 5.6** (contribution of far away assignments). Assume $n \to \infty$ and $d \to \infty$, and choose $t_{nd}$ decreasing to 0 such that $t_{nd}^2 \gg \frac{n+d}{nd}$. Then conditionally on $\Omega_1$ and for $n, d$ large enough that $2\sqrt{2nd}t_{nd} \ge gm$, we have:

$$ \sup_{\theta \in \Theta} \sum_{(\mathbf{z}, \mathbf{w}) \notin S(\mathbf{z}^*, \mathbf{w}^*, t_{nd})} p(\mathbf{z}, \mathbf{w}, \mathbf{x}; \theta) = p(\mathbf{z}^*, \mathbf{w}^*, \mathbf{x}; \theta^*) o_P(1) $$

where the $o_P$ is uniform in probability over all $\theta \in \Theta$.

## 5.2. Local Control

Proposition 5.5 gives deviations of order $O_P(\sqrt{nd})$, which are only useful for $(\mathbf{z}, \mathbf{w})$ such that $G$ and $\tilde{\Lambda}$ are large compared to $\sqrt{nd}$. For $(\mathbf{z}, \mathbf{w})$ close to $(\mathbf{z}^*, \mathbf{w}^*)$, we need tighter concentration inequalities, of order $o_P(n+d)$, as follows:

**Proposition 5.7** (small deviations $F_{nd}$). Conditionally upon $\Omega_1$, for $n$ and $d$ satisfying $(H_4)$, and for $(\mathbf{z}, \mathbf{w}) \in S(\mathbf{z}^*, \mathbf{w}^*, c/4)$ where $c > 0$ is defined in $(H_1)$, we have:

$$ \sum_{\substack{(\mathbf{z}, \mathbf{w}) \in S(\mathbf{z}^*, \mathbf{w}^*, c/4) \\ (\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)}} \frac{\Lambda(\mathbf{z}, \mathbf{w}) - \tilde{\Lambda}(\mathbf{z}^*, \mathbf{w}^*)}{d \|\mathbf{z} - \mathbf{z}^*\|_{0,\sim} + n \|\mathbf{w} - \mathbf{w}^*\|_{0,\sim}} = o_P(1) $$
The next proposition uses Propositions 4.6 and 5.7 to show that the combined contribution to the observed likelihood of assignments close to $(\mathbf{z}^*, \mathbf{w}^*)$ is also a $o_P$ of $p(\mathbf{z}^*, \mathbf{w}^*, \mathbf{x}; \mathbf{\theta}^*)$:

**Proposition 5.8 (contribution of local assignments).** With the previous notations and for $n$ and $d$ satisfying Assumption (H$_4$), for any $\tilde{c} \le \min(C, c/4)$ we have:

$$ \sup_{\mathbf{\theta} \in \Theta} \sum_{\substack{(\mathbf{z}, \mathbf{w}) \in S(\mathbf{z}^*, \mathbf{w}^*, \tilde{c}) \\ (\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)}} p(\mathbf{z}, \mathbf{w}, \mathbf{x}; \mathbf{\theta}) = o_P(p(\mathbf{z}^*, \mathbf{w}^*, \mathbf{x}; \mathbf{\theta}^*)) $$

### 5.3. Equivalent assignments

It remains to study the contribution of equivalent assignments.

**Proposition 5.9 (contribution of equivalent assignments).** For all $\mathbf{\theta} \in \Theta$, we have

$$ \sum_{(\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)} \frac{p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \mathbf{\theta})}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \mathbf{\theta}^*)} = \# \mathrm{Sym}(\mathbf{\theta}) \max_{\mathbf{\theta}' \sim \mathbf{\theta}} \frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \mathbf{\theta}')}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \mathbf{\theta}^*)} (1 + o_P(1)) $$

where the $o_P$ is uniform in $\mathbf{\theta}$.

The maximum over $\mathbf{\theta}' \sim \mathbf{\theta}$ accounts for equivalent configurations whereas $\# \mathrm{Sym}(\mathbf{\theta})$ is needed when $\mathbf{\theta}$ exhibits symmetry, as noticed in Remark 5.3.

## 6. Asymptotics for the Maximum Likelihood (MLE) and Variational (VE) Estimators

This section is devoted to the asymptotics of the MLE and VE in the incomplete data model as a consequence of the main result 5.1.

### 6.1. ML estimator

**Theorem 6.1 (Asymptotic behavior of $\hat{\theta}_{MLE}$).** Denote $\hat{\theta}_{MLE}$ the maximum likelihood estimator and use the notations of Proposition 3.1. There exist per-mutations $s$ of $\{1, \dots, g\}$ and $t$ of $\{1, \dots, m\}$ such that

$$ \hat{\pi}(\mathbf{z}^*) - \hat{\pi}_{MLE}^s = o_P(n^{-1/2}), \qquad \hat{\rho}(\mathbf{w}^*) - \hat{\rho}_{MLE}^t = o_P(d^{-1/2}), \\ \hat{\alpha}(\mathbf{z}^*, \mathbf{w}^*) - \hat{\alpha}_{MLE}^{s,t} = o_P((nd)^{-1/2}). $$

The proof relies on a Taylor expansion of the complete likelihood near its optimum, like in Proposition 3.2, and on our main theorem.
**Proof.**

Note first that unless $\Theta$ is constrained and with high probability, $\hat{\theta}_{MLE}$ and $\hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*)$ exhibit no symmetries. Indeed, equalities like $\hat{x}_{k\ell} = \hat{x}_{k',\ell'}$ have vanishingly small probabilities of being simultaneously true when $X_{ij}$ is discrete and null when $X_{ij}$ is continuous.

Note also that $\mathcal{L}_c^*$ has a unique maximum at $\hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*)$. Furthermore, the curvature of $\mathcal{L}_c^*$ at $\hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*)$ with respect to $\pi$ (resp. $\rho, \alpha$) converge in probability to $I_{\pi^*}$ (resp. $I_{\rho^*}, I_{\alpha^*}$) defined in Proposition 3.2 by consistency of $\hat{\theta}_c$. Therefore any estimator $\hat{\theta}$ bounded away from $\hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*)$ satisfies $\mathcal{L}_c^*(\hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*)) - \mathcal{L}_c^*(\hat{\theta}) > \Omega_P(1)$. If $\|\hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*) - \hat{\theta}\| = o_P(1)$, a Taylor expansion at $\hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*)$ gives

$$ \mathcal{L}_c^*(\hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*)) - \mathcal{L}_c^*(\hat{\theta}) = (n(\hat{\pi}_c(\mathbf{z}^*, \mathbf{w}^*) - \hat{\pi})^T I_{\pi^*} (\hat{\pi}_c(\mathbf{z}^*, \mathbf{w}^*) - \hat{\pi}) + \\ d(\hat{\rho}_c(\mathbf{z}^*, \mathbf{w}^*) - \hat{\rho})^T I_{\rho^*} (\hat{\rho}_c(\mathbf{z}^*, \mathbf{w}^*) - \hat{\rho}) + \\ nd \operatorname{Tr} [\{(\hat{\alpha}_c(\mathbf{z}^*, \mathbf{w}^*) - \hat{\alpha}) \odot (\hat{\alpha}_c(\mathbf{z}^*, \mathbf{w}^*) - \hat{\alpha})\}^T I_{\alpha^*}] \\ \times (1 + o_P(1)) + o_P(1) $$

where the linear term in the expansion vanishes as $\hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*)$ is the argmax of $\mathcal{L}_c^*$ and the Hessian of $\mathcal{L}_c^*$ at $\hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*)$ were replaced by their limit in probability.

We may now prove the corollary by contradiction. Assume that $\min_s(\hat{\pi}_{MLE}^s - \hat{\pi}(\mathbf{w}^*)) \neq o_P(n^{-1/2})$, $\min_t(\hat{\rho}_{MLE}^t - \hat{\rho}(\mathbf{w}^*)) \neq o_P(d^{-1/2})$ or $\min_{s,t}(\hat{\alpha}_{MLE}^{s,t} - \hat{\alpha}(\mathbf{z}^*, \mathbf{w}^*)) \neq o_P(n d^{-1/2})$ where $s$ and $t$ are permutations of $\{1, \dots, g\}$ and $\{1, \dots, m\}$. Plugging $\hat{\theta}_{MLE}^{s,t}$ in the previous expansion shows that:

$$ \min_{s,t} \mathcal{L}_c^*(\hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*)) - \mathcal{L}_c^*(\hat{\theta}_{MLE}^{s,t}) = \Omega_P(1). \quad (6.1) $$

But, since $\hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*)$ and $\hat{\theta}_{MLE}$ maximize respectively $\frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta)}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*)}$ and $\frac{p(\mathbf{x}; \theta)}{p(\mathbf{x}; \theta^*)}$ and have no symmetries, it follows by Theorem 5.1 that

$$ \left| \frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \hat{\theta}_c(\mathbf{z}^*, \mathbf{w}^*))}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*)} - \max_{s,t} \frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \hat{\theta}_{MLE}^{s,t})}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*)} \right| = o_P(1) $$

which contradicts Equation (6.1) and concludes the proof.

□

## 6.2. Variational estimator

Due to the complex dependence structure of the observations, the maximum likelihood estimator of the LBM is not numerically tractable, even with the EM-algorithm. In practice, a variational approximation can be used, see for
example [5]: for any joint distribution $Q \in Q$ on $Z \times W$ a lower bound of $\mathcal{L}(\theta)$ is given by

$$ 
\begin{aligned}
J(Q, \theta) &= \mathcal{L}(\theta) - KL(Q, p(\cdot, \cdot; \theta, \mathbf{x})) \\
&= \mathbb{E}_Q[\mathcal{L}_c(\mathbf{z}, \mathbf{w}; \theta)] + \mathcal{H}(Q).
\end{aligned}
$$

where $\mathcal{H}(Q) = -\mathbb{E}_Q[\log(Q)]$. Choose $Q$ to be the set of factorized distributions, such that for all $(\mathbf{z}, \mathbf{w})$

$$ \mathbb{Q}(\mathbf{z}, \mathbf{w}) = \mathbb{Q}(\mathbf{z}) \mathbb{Q}(\mathbf{w}) = \prod_{i,k} \mathbb{Q}(z_{ik} = 1)^{z_{ik}} \prod_{j,\ell} \mathbb{Q}(w_{j\ell} = 1)^{w_{j\ell}} $$

allows to obtain tractable expressions of $J(Q, \theta)$ as a lower bound of the log-likelihood. The variational estimate $\hat{\theta}_{var}$ of $\theta$ is defined as

$$ \hat{\theta}_{var} \in \underset{\theta \in \Theta}{\operatorname{argmax}} \max_{Q \in \mathbb{Q}} J(Q, \theta). $$

The following corollary states that $\hat{\theta}_{var}$ has the same asymptotics as $\hat{\theta}_{MLE}$ and $\hat{\theta}_{c}$.

**Theorem 6.2 (Variational estimate).** Under the assumptions of Theorem 5.1 there exist permutations $s$ of $\{1, \dots, g\}$ and $t$ of $\{1, \dots, m\}$ such that

$$ 
\begin{aligned}
\hat{\pi}(\mathbf{z}^*) - \hat{\pi}_{var}^s &= o_P(n^{-1/2}), & \hat{\rho}(\mathbf{w}^*) - \hat{\rho}_{var}^t &= o_P(d^{-1/2}), \\
\hat{\alpha}(\mathbf{z}^*, \mathbf{w}^*) - \hat{\alpha}_{var}^{s,t} &= o_P((nd)^{-1/2}).
\end{aligned}
$$

**Proof.**

Remark first that for every $\theta$ and for every $(\mathbf{z}, \mathbf{w})$,

$$ p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \theta) \le \exp[J(\delta_{\mathbf{z}} \times \delta_{\mathbf{w}}, \theta)] \le \max_{Q \in \mathbb{Q}} \exp[J(Q, \theta)] \le p(\mathbf{x}; \theta) $$

where $\delta_{\mathbf{z}}$ denotes the dirac mass on $\mathbf{z}$. By dividing by $p(\mathbf{x}; \theta^*)$, we obtain

$$ \frac{p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \theta)}{p(\mathbf{x}; \theta^*)} \le \frac{\max_{Q \in \mathbb{Q}} \exp[J(Q, \theta)]}{p(\mathbf{x}; \theta^*)} \le \frac{p(\mathbf{x}; \theta)}{p(\mathbf{x}; \theta^*)}. $$

As this inequality is true for every couple $(\mathbf{z}, \mathbf{w})$, we have in particular:

$$ \max_{(\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)} \frac{p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \theta)}{p(\mathbf{x}; \theta^*)} = \max_{\theta' \sim \theta} \frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta')}{p(\mathbf{x}; \theta^*)} \le \frac{\max_{Q \in \mathbb{Q}} \exp[J(Q, \theta)]}{p(\mathbf{x}; \theta^*)}. $$

Noticing that $p(\mathbf{x}; \theta^*) = \#\mathrm{Sym}(\theta^*)p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*) (1 + o_p(1))$, Theorem 5.1 therefore leads to the following bounds:

$$ 
\begin{aligned}
\max_{\theta' \sim \theta} \frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta')}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*)} (1 + o_P(1)) &\le \frac{\max_{Q \in \mathbb{Q}} \exp[J(Q, \theta)]}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*)} \\
&\le \#\mathrm{Sym}(\theta) \max_{\theta' \sim \theta} \frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta')}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*)} (1 + o_P(1)) + o_P(1).
\end{aligned}
$$
Again, unless $\Theta$ is constrained, $\hat{\theta}_{VAR}$ exhibits no symmetries with high probability and the same proof by contradiction as in section 6.1 gives the result. $\square$

## 7. Conclusion

The Latent Block Model offers challenging theoretical questions. We solved under mild assumptions the consistency and asymptotic normality of the maximum likelihood and variational estimators for observations with conditional density belonging to a univariate exponential family, and for a balanced asymptotic rate between the number of rows $n$ and the number of columns $d$: $\log(d)/n \to 0$ and $\log(n)/d \to 0$ as $n$ and $d$ tend to infinity. Our results extend those of [2] for binary SBM not only by managing the double direction of LBM, but also by considering larger types of observations. That brought us to define specific concentration inequalities as large and moderate deviations concerning sub-exponential variables. Moreover, we dealt with specific cases of symmetry that were not taken into account as of now.

A specific framework of sparsity was studied by [2]. This is especially convenient for SBM, to model reasonable network settings: increasing sparsity (i.e. number of 0) can be done directly by scaling the Bernoulli parameters $p_{k\ell}$ with a common factor that should decrease no faster than $\Omega(\log^\delta(n)/n)$, with $\delta > 2$, to ensure consistency. This could also be considered for binary LBM. However this approach fails to model actual observations in the more general valued setting. The equivalent approach could be to consider the product of a Bernoulli variable with the actual observation value. Note however, than even without considering sparsity we recover essentially the same rate: in the sparse-SBM case, each node should be connected to $\Omega(\log^\delta(n))$ others to ensure consistency whereas in the dense-LBM case, each of the $n$-row should should be characterized by $\Omega(\log^\delta(n))$ columns (and vice-versa) to ensure consistency.

Alternative research direction could be to explore asymptotic settings where the numbers $n$ of rows and $d$ columns grow at very different rates. Other open question concern estimation of the number of row and column groups and settings where the number of groups increases with $n$ and $d$.
## Appendix A: Proofs of section 4

### A.1. Proof of Proposition 4.5 (maximum of G and $\tilde{\Lambda}$ in $\theta$)

**Proof.**

Define $\nu(x, \alpha) = x\alpha - \psi(\alpha)$. For $x$ fixed, $\nu(x, \alpha)$ is maximized at $\alpha = (\psi')^{-1}(x)$. Manipulations yield

$$F_{nd}(\alpha, \mathbf{z}, \mathbf{w}) = \log p(\mathbf{x}; \mathbf{z}, \mathbf{w}, \theta) - \log p(\mathbf{x}; \mathbf{z}^*, \mathbf{w}^*, \theta^*) \\ = nd \left[ \sum_k \sum_l \hat{\pi}_k(\mathbf{z}) \hat{\rho}_l(\mathbf{w}) \nu(\hat{x}_{k\ell}(\mathbf{z}, \mathbf{w}), \alpha_{k\ell}) - \sum_k \sum_l \hat{\pi}_k(\mathbf{z}^*) \hat{\rho}_l(\mathbf{w}^*) \nu(\hat{x}_{k\ell}(\mathbf{z}^*, \mathbf{w}^*), \alpha_{k\ell}^*) \right]$$

which is maximized at $\alpha_{k\ell} = (\psi')^{-1}(\hat{x}_{k\ell}(\mathbf{z}, \mathbf{w}))$. Similarly

$$G(\alpha, \mathbf{z}, \mathbf{w}) = \mathbb{E}_{\theta^*} \left[ \log p(\mathbf{x}; \mathbf{z}, \mathbf{w}, \theta) - \log p(\mathbf{x}; \mathbf{z}^*, \mathbf{w}^*, \theta^*) | \mathbf{z}^*, \mathbf{w}^* \right] \\ = nd \left[ \sum_k \sum_l \hat{\pi}_k(\mathbf{z}) \hat{\rho}_l(\mathbf{w}) \nu(\bar{x}_{k\ell}(\mathbf{z}, \mathbf{w}), \alpha_{k\ell}) - \sum_k \sum_l \hat{\pi}_k(\mathbf{z}^*) \hat{\rho}_l(\mathbf{w}^*) \nu(\psi'(\alpha_{k\ell}^*), \alpha_{k\ell}^*) \right]$$

is maximized at $\alpha_{k\ell} = (\psi')^{-1}(\bar{x}_{k\ell}(\mathbf{z}, \mathbf{w}))$

<figure><img src="image_6.png" /></figure>

### A.2. Proof of Proposition 4.6 (maximum of G and $\tilde{\Lambda}$ in $(\theta, \mathbf{z}, \mathbf{w})$)

**Proof.**

We condition on $(\mathbf{z}^*, \mathbf{w}^*)$ and prove Equation (4.2):

$$G(\theta, \mathbf{z}, \mathbf{w}) = \mathbb{E}_{\theta^*} \left[ \log \frac{p(\mathbf{x}; \mathbf{z}, \mathbf{w}, \theta)}{p(\mathbf{x}; \mathbf{z}^*, \mathbf{w}^*, \theta^*)} \middle| \mathbf{z}^*, \mathbf{w}^* \right] \\ = \sum_i \sum_j \sum_{k,k'} \sum_{\ell, \ell'} \mathbb{E}_{\theta^*} [x_{ij}(\alpha_{k'\ell'} - \alpha_{k\ell}^*) - (\psi(\alpha_{k'\ell'}) - \psi(\alpha_{k\ell}^*))] z_{ik}^* z_{ik'} w_{j\ell}^* w_{j\ell'} \\ = nd \sum_{k,k'} \sum_{\ell, \ell'} \mathbb{R}_g(\mathbf{z})_{k,k'} \mathbb{R}_m(\mathbf{w})_{\ell,\ell'} [\psi'(\alpha_{k\ell}^*)(\alpha_{k'\ell'} - \alpha_{k\ell}^*) + \psi(\alpha_{k\ell}^*) - \psi(\alpha_{k'\ell'})] \\ = -nd \sum_{k,k'} \sum_{\ell, \ell'} \mathbb{R}_g(\mathbf{z})_{k,k'} \mathbb{R}_m(\mathbf{w})_{\ell,\ell'} \text{KL}(\alpha_{k\ell}^*, \alpha_{k'\ell'})$$

If $(\mathbf{z}^*, \mathbf{w}^*)$ is regular, and for $n, d > 2/c$, all the rows of $\mathbb{R}_g(\mathbf{z})$ and $\mathbb{R}_m(\mathbf{w})$ have at least one positive element and we can apply lemma D.1 (which is an adaptation for LBM of Lemma 3.2 of [2] for SBM) to characterize the maximum for G.

The maximality of $\tilde{\Lambda}(\mathbf{z}^*, \mathbf{w}^*)$ results from the fact that $\tilde{\Lambda}(\mathbf{z}, \mathbf{w}) = G(\tilde{\alpha}(\mathbf{z}, \mathbf{w}), \mathbf{z}, \mathbf{w})$ where $\tilde{\alpha}(\mathbf{z}, \mathbf{w})$ is a particular value of $\alpha$, $\tilde{\Lambda}$ is immediately maximum at $(\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)$, and for those, we have $\tilde{\alpha}(\mathbf{z}, \mathbf{w}) \sim \alpha^*$.

The separation and local behavior of G around $(\mathbf{z}^*, \mathbf{w}^*)$ is a direct consequence of the proposition 4.7.

<figure><img src="image_7.png" /></figure>
### A.3. Proof of Proposition 4.7 (Local upper bound for $\tilde{\Lambda}$)

**Proof.**

We work conditionally on $(\mathbf{z}^*, \mathbf{w}^*)$. The principle of the proof relies on the extension of $\tilde{\Lambda}$ to a continuous subspace of $M_g([0, 1]) \times M_m([0, 1])$, in which confusion matrices are naturally embedded. The regularity assumption allows us to work on a subspace that is bounded away from the borders of $M_g([0, 1]) \times M_m([0, 1])$. The proof then proceeds by (1) computing the gradient of $\tilde{\Lambda}$ at and around its argmax and (2) using those gradients to control the local behavior of $\tilde{\Lambda}$ around its argmax. The local behavior allows in turn to show that $\tilde{\Lambda}$ is well-separated.

Note that $\tilde{\Lambda}$ only depends on $\mathbf{z}$ and $\mathbf{w}$ through $\mathbb{R}_g(\mathbf{z})$ and $\mathbb{R}_m(\mathbf{w})$. We can therefore extend it to matrices $(U, V) \in \mathcal{U}_c \times \mathcal{V}_c$ where $\mathcal{U}$ is the subset of matrices $M_g([0, 1])$ with each row sum higher than $c/2$ and $\mathcal{V}$ is a similar subset of $M_m([0, 1])$.

$$ \tilde{\Lambda}(U, V) = -nd \sum_{k,k'} \sum_{\ell, \ell'} U_{kk'} V_{\ell\ell'} \text{KL}(\alpha_{k\ell}^*, \bar{\alpha}_{k'\ell'}) $$

where

$$ \bar{\alpha}_{k\ell} = \bar{\alpha}_{k\ell}(U, V) = (\psi')^{-1} \left( \frac{[U^T S^* V]_{k\ell}}{[U^T 1 V]_{k\ell}} \right) $$

and $1$ is the $g \times m$ matrix filled with 1. Confusion matrices $\mathbb{R}_g(\mathbf{z})$ and $\mathbb{R}_m(\mathbf{w})$ satisfy $\mathbb{R}_g(\mathbf{z})\mathbb{I} = \hat{\pi}(\mathbf{z}^*)$ and $\mathbb{R}_m(\mathbf{w})\mathbb{I} = \hat{\rho}(\mathbf{w}^*)$, with $\mathbb{I} = (1, \dots, 1)^T$ a vector only containing 1 values, and are obviously in $\mathcal{U}_c$ and $\mathcal{V}_c$ as soon as $(\mathbf{z}^*, \mathbf{w}^*)$ is $c/2$ regular.

The maps $f_{k,\ell} : (U, V) \mapsto KL(\alpha_{k\ell}^*, \bar{\alpha}_{k\ell}(U, V))$ are twice differentiable with second derivatives bounded over $\mathcal{U}_c \times \mathcal{V}_c$ and therefore so is $\tilde{\Lambda}(U, V)$. Tedious but straightforward computations show that the derivative of $\tilde{\Lambda}$ at $(D_\pi, D_\rho) := (\text{Diag}(\hat{\pi}(\mathbf{z}^*)), \text{Diag}(\hat{\rho}(\mathbf{w}^*)))$ is:

$$ A_{kk'}(\mathbf{w}^*) := -\frac{1}{nd} \frac{\partial \tilde{\Lambda}}{\partial U_{kk'}}(D_\pi, D_\rho) = \sum_{\ell} \hat{\rho}_\ell(\mathbf{w}^*) \text{KL}(\alpha_{k\ell}^*, \alpha_{k'\ell}^*) $$

$$ B_{\ell\ell'}(\mathbf{z}^*) := -\frac{1}{nd} \frac{\partial \tilde{\Lambda}}{\partial V_{\ell\ell'}}(D_\pi, D_\rho) = \sum_k \widehat{\pi}_k(\mathbf{z}^*) \text{KL}(\alpha_{k\ell}^*, \alpha_{k\ell'}^*) $$

$A(\mathbf{w}^*)$ and $B(\mathbf{z}^*)$ are the matrix-derivative of $-\tilde{\Lambda}/nd$ at $(D_\pi, D_\rho)$. Since $(\mathbf{z}^*, \mathbf{w}^*)$ is $c/2$-regular and by definition of $\delta(\alpha^*)$, $A(\mathbf{w}^*)_{kk'} \ge c\delta(\alpha^*)/2$ (resp. $B(\mathbf{w}^*)_{\ell\ell'} \ge c\delta(\alpha^*)/2$) if $k \ne k'$ (resp. $\ell \ne \ell'$) and $A(\mathbf{w}^*)_{kk} = 0$ (resp. $B(\mathbf{z}^*)_{\ell\ell} = 0$) for all $k$ (resp. $\ell$). By boundedness of the second derivative, there exists $C > 0$ such that for all $(D_\pi, D_\rho)$ and all $(H, G) \in S(D_\pi, D_\rho, C)$, where the definition of the set of local assignments $S$ is extended to the subset of matrices, we have:

$$ -\frac{1}{nd} \frac{\partial \tilde{\Lambda}}{\partial U_{kk'}}(H, G) \left\{ \begin{array}{ll} \ge \frac{3c\delta(\alpha^*)}{8} & \text{if } k \ne k' \\ \le \frac{c\delta(\alpha^*)}{8} & \text{if } k = k' \end{array} \right. \quad \text{and} \quad -\frac{1}{nd} \frac{\partial \tilde{\Lambda}}{\partial V_{\ell\ell'}}(H, G) \left\{ \begin{array}{ll} \ge \frac{3c\delta(\alpha^*)}{8} & \text{if } \ell \ne \ell' \\ \le \frac{c\delta(\alpha^*)}{8} & \text{if } \ell = \ell' \end{array} \right. $$
Choose $U$ and $V$ in $(\mathcal{U}_c \times \mathcal{V}_c) \cap S(D_\pi, D_\rho, C)$ satisfying $U\mathbb{I} = \pi(\mathbf{z}^*)$ and $V\mathbb{I} = \rho(\mathbf{w}^*)$. $U - D_\pi$ and $V - D_\rho$ have nonnegative off diagonal coefficients and negative diagonal coefficients. Furthermore, the coefficients of $U, V, D_\pi, D_\rho$ sum up to 1 and $\mathrm{Tr}(D_\pi) = \mathrm{Tr}(D_\rho) = 1$. By Taylor expansion, there exists a couple $(H, G)$ also in $(\mathcal{U}_c \times \mathcal{V}_c) \cap S(D_\pi, D_\rho, C)$ such that

$$ \frac{-1}{nd} \tilde{\Lambda}(U, V) = \frac{-1}{nd} \left[ \tilde{\Lambda}(D_\pi, D_\rho) + \mathrm{Tr} \left( (U - D_\pi) \frac{\partial \tilde{\Lambda}}{\partial U}(H, G) \right) + \mathrm{Tr} \left( (V - D_\rho) \frac{\partial \tilde{\Lambda}}{\partial V}(H, G) \right) \right] \\ \geq \frac{c\delta(\alpha^*)}{8} \left[ 3 \sum_{k \neq k'} (U - D_\pi)_{kk'} + 3 \sum_{\ell \neq \ell'} (V - D_\rho)_{\ell\ell'} + \sum_k (U - D_\pi)_{kk} + \sum_\ell (V - D_\rho)_{\ell\ell} \right] \\ = \frac{c\delta(\alpha^*)}{4} [(1 - \mathrm{Tr}(U)) + (1 - \mathrm{Tr}(V))] $$

To conclude the proof, assume without loss of generality that $(\mathbf{z}, \mathbf{w}) \in S(\mathbf{z}^*, \mathbf{w}^*, C)$ achieves the $\|\cdot\|_{0,\sim}$ norm (i.e. it is the closest to $(\mathbf{z}^*, \mathbf{w}^*)$ in its representative class). Then $(U, V) = (\mathbb{R}_g(\mathbf{z}), \mathbb{R}_m(\mathbf{w}))$ is in $(\mathcal{U}_c \times \mathcal{V}_c) \cap S(D_\pi, D_\rho, C)$ and satisfy $U\mathbb{I} = \pi(\mathbf{z}^*)$ (resp. $V\mathbb{I} = \rho(\mathbf{w}^*)$). We just need to note $n(1 - \mathrm{Tr}(\mathbb{R}_g(\mathbf{z}))) = \|\mathbf{z} - \mathbf{z}^*\|_{0,\sim}$ (resp. $d(1 - \mathrm{Tr}(\mathbb{R}_m(\mathbf{w}))) = \|\mathbf{w} - \mathbf{w}^*\|_{0,\sim}$) to end the proof. $\square$

<figure><img src="image_8.png" /></figure>

# Appendix B: Proofs of section 5

## B.1. Proof of Proposition 5.5 (large deviation for $F_{nd}$)

### Proof.

Conditionally upon $(\mathbf{z}^*, \mathbf{w}^*)$,

$$ F_{nd}(\theta, \mathbf{z}, \mathbf{w}) - \tilde{\Lambda}(\mathbf{z}, \mathbf{w}) \leq F_{nd}(\theta, \mathbf{z}, \mathbf{w}) - G(\theta, \mathbf{z}, \mathbf{w}) \\ = \sum_i \sum_j (\alpha_{z_i w_j} - \alpha_{z_i^* w_j^*}^*) (x_{ij} - \psi'(\alpha_{z_i^* w_j^*}^*)) \\ = \sum_{kk'} \sum_{\ell\ell'} (\alpha_{k'\ell'} - \alpha_{k\ell}^*) W_{kk'\ell\ell'} \\ \leq \sup_{\Gamma \in \mathbb{R}^{g^2 \times m^2} \atop \|\Gamma\|_{\infty} \leq \mathrm{Diam}(\Theta)} \sum_{kk'} \sum_{\ell\ell'} \Gamma_{kk'\ell\ell'} W_{kk'\ell\ell'} := Z $$

uniformly in $\theta$, where the $W_{kk'\ell\ell'}$ are independent and defined by:

$$ W_{kk'\ell\ell'} = \sum_i \sum_j z_{ik}^* w_{j\ell}^* z_{i,k'} w_{j\ell'} (x_{ij} - \psi'(\alpha_{k\ell}^*)) $$
is the sum of $nd\mathbb{R}_g(\mathbf{z})_{kk'}\mathbb{R}_m(\mathbf{w})_{\ell\ell'}$ sub-exponential variables with parameters $(\bar{\sigma}^2, 1/\kappa)$ and is therefore itself sub-exponential with parameters $(nd\mathbb{R}_g(\mathbf{z})_{kk'}\mathbb{R}_m(\mathbf{w})_{\ell\ell'}\bar{\sigma}^2, 1/\kappa)$. According to Proposition C.4, $\mathbb{E}_{\mathbf{\theta}^*}[Z|\mathbf{z}^*, \mathbf{w}^*] \le gm \text{Diam}(\mathbf{\Theta})\sqrt{nd\bar{\sigma}^2}$ and $Z$ is sub-exponential with parameters $(nd \text{Diam}(\mathbf{\Theta}))^2(2\sqrt{2})^2\bar{\sigma}^2, 2\sqrt{2} \text{Diam}(\mathbf{\Theta})/\kappa)$. In particular, for all $\varepsilon_{n,d} < \bar{\sigma}\kappa$.

$$ \begin{aligned} \mathbb{P}_{\mathbf{\theta}^*} \left( Z \ge \bar{\sigma}gm \text{Diam}(\mathbf{\Theta})\sqrt{nd} \left\{ 1 + \frac{\sqrt{8nd\varepsilon_{n,d}}}{gm} \right\} \middle| \mathbf{z}^*, \mathbf{w}^* \right) & \\ \le \mathbb{P}_{\mathbf{\theta}^*} \left( Z \ge \mathbb{E}_{\mathbf{\theta}^*}[Z|\mathbf{z}^*, \mathbf{w}^*] + \bar{\sigma} \text{Diam}(\mathbf{\Theta})nd2\sqrt{2}\varepsilon_{n,d} \middle| \mathbf{z}^*, \mathbf{w}^* \right) & \\ & \le \exp \left( - \frac{nd\varepsilon_{n,d}^2}{2} \right) \end{aligned} $$

We can then remove the conditioning and take a union bound to prove Equation (5.1). □

## B.2. Proof of Proposition 5.6 (contribution of far away assignments)

### Proof.

Conditionally on $(\mathbf{z}^*, \mathbf{w}^*)$, we know from proposition 4.6 that $\tilde{\Lambda}$ is maximal in $(\mathbf{z}^*, \mathbf{w}^*)$ and its equivalence class. Choose $0 < t_{nd}$ decreasing to 0 but satisfying $t_{nd}^2 \gg \frac{n+d}{nd}$. This is possible as $n \to 0$ and $d \to 0$.

According to equation 4.4, for all $(\mathbf{z}, \mathbf{w}) \notin S(\mathbf{z}^*, \mathbf{w}^*, C)$

$$ \tilde{\Lambda}(\mathbf{z}, \mathbf{w}) \le -B(C) \text{nd.} $$

Now, according to Equation 4.3, for all $(\mathbf{z}, \mathbf{w}) \in S(\mathbf{z}^*, \mathbf{w}^*, C) \setminus S(\mathbf{z}^*, \mathbf{w}^*, t_{nd})$

$$ \tilde{\Lambda}(\mathbf{z}, \mathbf{w}) \le -\frac{c\delta(\boldsymbol{\alpha}^*)}{4} (n\|\mathbf{w} - \mathbf{w}^*\|_{0,\sim} + d\|\mathbf{z} - \mathbf{z}^*\|_{0,\sim}) \le -\frac{c\delta(\boldsymbol{\alpha}^*)}{4} ndt_{nd} $$

since either $\|\mathbf{z} - \mathbf{z}^*\|_{0,\sim} \ge nt_{nd}$ or $\|\mathbf{w} - \mathbf{w}^*\|_{0,\sim} \ge dt_{nd}$. Hence, for $n$ and $d$ large enough for all $(\mathbf{z}, \mathbf{w}) \notin S(\mathbf{z}^*, \mathbf{w}^*, t_{nd})$

$$ \tilde{\Lambda}(\mathbf{z}, \mathbf{w}) \le -\frac{c\delta(\boldsymbol{\alpha}^*)}{4} ndt_{nd} \quad (B.1) $$

Set $\varepsilon_{nd} = \inf \left( \frac{c\delta(\boldsymbol{\alpha}^*)t_{nd}}{32\sqrt{2\bar{\sigma}} \text{Diam}(\mathbf{\Theta})}, \kappa\bar{\sigma} \right)$. By proposition 5.5, and with our choice of $\varepsilon_{nd}$, with probability higher than $1 - \Delta_{nd}^1(\varepsilon_{nd})$,
$$
\begin{align*}
& \sum_{(\mathbf{z}, \mathbf{w}) \notin S(\mathbf{z}^*, \mathbf{w}^*, t_{nd})} p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \boldsymbol{\theta}) \\
&= p(\mathbf{x}|\mathbf{z}^*, \mathbf{w}^*, \boldsymbol{\theta}^*) \sum_{(\mathbf{z}, \mathbf{w}) \notin S(\mathbf{z}^*, \mathbf{w}^*, t_{nd})} p(\mathbf{z}, \mathbf{w}; \boldsymbol{\theta}) e^{F_{nd}(\boldsymbol{\theta}, \mathbf{z}, \mathbf{w}) - \tilde{\Lambda}(\mathbf{z}, \mathbf{w}) + \tilde{\Lambda}(\mathbf{z}, \mathbf{w})} \\
&\le p(\mathbf{x}|\mathbf{z}^*, \mathbf{w}^*, \boldsymbol{\theta}^*) \sum_{(\mathbf{z}, \mathbf{w}) \in \mathcal{Z} \times \mathcal{W}} p(\mathbf{z}, \mathbf{w}; \boldsymbol{\theta}) e^{F_{nd}(\boldsymbol{\theta}, \mathbf{z}, \mathbf{w}) - \tilde{\Lambda}(\mathbf{z}, \mathbf{w}) - ndt_{nd}c\delta(\boldsymbol{\alpha}^*)/4} \\
&\le p(\mathbf{x}|\mathbf{z}^*, \mathbf{w}^*, \boldsymbol{\theta}^*) \sum_{\mathbf{z}, \mathbf{w}} p(\mathbf{z}, \mathbf{w}; \boldsymbol{\theta}) e^{-ndt_{nd}c\delta(\boldsymbol{\alpha}^*)/8} \\
&= \frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*)}{p(\mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*)} e^{-ndt_{nd}c\delta(\boldsymbol{\alpha}^*)/8} \\
&\le p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*) \exp \left( -ndt_{nd} \frac{c\delta(\boldsymbol{\alpha}^*)}{8} + (n+d) \log \frac{1}{c} \right) \\
&= p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*) o(1)
\end{align*}
$$

where the second line comes from inequality (B.1), the third from the global control studied in Proposition 5.5 and the definition of $\varepsilon_{nd}$, the fourth from the definition of $p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*)$, the fifth from the bounds on $\pi^*$ and $\rho^*$ and the last from $t_{nd} \gg (n+d)/nd$.

In addition, we have $\varepsilon_{nd}^2 \gg \frac{n+d}{nd}$ so that $\Delta_{nd}^1(\varepsilon_{nd})$ vanishes and:

$$ \sum_{(\mathbf{z}, \mathbf{w}) \notin S(\mathbf{z}^*, \mathbf{w}^*, t_{nd})} p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \boldsymbol{\theta}) = p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*, \boldsymbol{\theta}^*) o_P(1) $$

where the $o_P$ is uniform in probability over all $\boldsymbol{\theta} \in \mathbf{\Theta}$.

<figure><img src="image_9.png" /></figure>

## B.3. Proof of Proposition 5.7 (local convergence $F_{nd}$)

### Proof.

We work conditionally on $(\mathbf{z}^*, \mathbf{w}^*) \in \mathcal{Z}_1 \times \mathcal{W}_1$. We assume with no loss of generality that $(\mathbf{z}, \mathbf{w})$ is the representative of its class closest to $(\mathbf{z}^*, \mathbf{w}^*)$. Choose $\varepsilon \le \kappa\underline{\sigma}^2$ small. Manipulation of $\tilde{\Lambda}$ and $\tilde{\Lambda}$ yield

$$
\begin{align*}
\frac{F_{nd}(\boldsymbol{\theta}, \mathbf{z}, \mathbf{w}) - \tilde{\Lambda}(\mathbf{z}, \mathbf{w})}{nd} &\le \frac{\Lambda(\mathbf{z}, \mathbf{w}) - \tilde{\Lambda}(\mathbf{z}, \mathbf{w})}{nd} \\
&= \sum_k \sum_\ell \hat{\pi}_k(\mathbf{z}) \hat{\rho}_\ell(\mathbf{w}) [f(\hat{x}_{k\ell}) - f(\bar{x}_{k\ell})] \\
&\quad - \sum_k \sum_\ell \hat{\pi}_k(\mathbf{z}^*) \hat{\rho}_\ell(\mathbf{w}^*) \alpha_{k\ell}^* (\hat{x}_{k\ell}^* - \bar{x}_{k\ell}^*)
\end{align*}
$$
where $f(x) = x(\psi')^{-1}(x) - \psi \circ (\psi')^{-1}(x)$, $\hat{x}_{k\ell}^* = \hat{x}_{k\ell}(\mathbf{z}^*, \mathbf{w}^*)$ and $\bar{x}_{k\ell}^* = \psi'(\alpha_{k\ell}^*)$. The function $f$ is twice differentiable on $\mathcal{A}$ with $f'(x) = -(\psi')^{-1}(x)$ and $f''(x) = -1/\psi'' \circ (\psi')^{-1}(x)$. $f'$ (resp. $f''$) are bounded over $I = \psi'([-M_\alpha, M_\alpha])$ by $M_\alpha$ (resp. $1/\sigma^2$).

Assignments $(\mathbf{z}, \mathbf{w})$ belonging to $S(\mathbf{z}^*, \mathbf{w}^*, c/4)$ are also $c/4$-regular. According to Proposition C.2, $\hat{x}_{k\ell}$ and $\bar{x}_{k\ell}$ are at distance at most $\varepsilon$ with probability higher than $1 - 2 \exp(-\frac{ndc^2\varepsilon^2}{32(\bar{\sigma}^2+\kappa^{-1}\varepsilon)})$, so that:

$$f(\hat{x}_{k\ell}) - f(\bar{x}_{k\ell}) = f'(\bar{x}_{k\ell})(\hat{x}_{k\ell} - \bar{x}_{k\ell}) + \Omega((\hat{x}_{k\ell} - \bar{x}_{k\ell})^2)$$

By Proposition C.2, $(\hat{x}_{k\ell} - \bar{x}_{k\ell})^2 = O_P(1/nd)$ where the $O_P$ is uniform in $\mathbf{z}, \mathbf{w}$ and does not depend on $\mathbf{z}^*, \mathbf{w}^*$. Similarly,

$$f'(\bar{x}_{k\ell}) = f'(\bar{x}_{k\ell}^*) + \Omega(\bar{x}_{k\ell} - \bar{x}_{k\ell}^*) = \alpha_{k\ell}^* + \Omega(\bar{x}_{k\ell} - \bar{x}_{k\ell}^*)$$

$\bar{x}_{k\ell}$ is a convex combination of the $S_{k\ell}^* = \psi'(\alpha_{k\ell}^*)$ therefore,

$$|\bar{x}_{k\ell} - \bar{x}_{k\ell}^*| = \left| \frac{[\mathbb{R}_g(\mathbf{z})^T \mathbf{S}^* \mathbb{R}_m(\mathbf{w})]_{k\ell}}{\hat{\pi}_k(\mathbf{z}) \hat{\rho}_\ell(\mathbf{w})} - \bar{x}_{k\ell}^* \right| \\ \leq \left( 1 - \frac{[\mathbb{R}_g(\mathbf{z})_{kk} \mathbb{R}_m(\mathbf{w})_{\ell\ell}]}{\hat{\pi}_k(\mathbf{z}) \hat{\rho}_\ell(\mathbf{w})} \right) (S_{\max}^* - S_{\min}^*)$$

Note that:

$$\sum_{k,\ell} \hat{\pi}_k(\mathbf{z}) \hat{\rho}_\ell(\mathbf{w}) \left( 1 - \frac{[\mathbb{R}_g(\mathbf{z})_{kk} \mathbb{R}_m(\mathbf{w})_{\ell\ell}]}{\hat{\pi}_k(\mathbf{z}) \hat{\rho}_\ell(\mathbf{w})} \right) = 1 - \mathrm{Tr}(\mathbb{R}_g(\mathbf{z})) \mathrm{Tr}(\mathbb{R}_m(\mathbf{w})) \\ \leq \frac{\|\mathbf{z} - \mathbf{z}^*\|_{0,\sim}}{n} + \frac{\|\mathbf{w} - \mathbf{w}^*\|_{0,\sim}}{d}$$

and $\hat{x}_{k\ell} - \bar{x}_{k\ell} = o_P(1)$. Therefore

$$\sum_{k,\ell} \hat{\pi}_k(\mathbf{z}) \hat{\rho}_\ell(\mathbf{w}) \Omega(\bar{x}_{k\ell} - \bar{x}_{k\ell}^*) \times (\hat{x}_{k\ell} - \bar{x}_{k\ell}) = o_P \left( \frac{\|\mathbf{z} - \mathbf{z}^*\|_{0,\sim}}{n} + \frac{\|\mathbf{w} - \mathbf{w}^*\|_{0,\sim}}{d} \right)$$

The remaining term writes

$$\sum_{k,\ell} \alpha_{k\ell}^* [\hat{\pi}_k(\mathbf{z}) \hat{\rho}_\ell(\mathbf{w})(\hat{x}_{k\ell} - \bar{x}_{k\ell}) - \hat{\pi}_k(\mathbf{z}^*) \hat{\rho}_\ell(\mathbf{w}^*)(\hat{x}_{k\ell}^* - \bar{x}_{k\ell}^*)]$$

According to Proposition C.3, this term is $o_P(\frac{\|\mathbf{z}-\mathbf{z}^*\|_{0,\sim}}{n} + \frac{\|\mathbf{w}-\mathbf{w}^*\|_{0,\sim}}{d})$ uniformly in $(\mathbf{z}, \mathbf{w})$ and $(\mathbf{z}^*, \mathbf{w}^*) \in \Omega_1$ as soon $ne^{-ad} \to 0$ and $de^{-an} \to 0$ for all $a > 0$, which is true under $(H_4)$. It follows that:

$$\sup_{\substack{(\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*) \\ (\mathbf{z}, \mathbf{w}) \in S(\mathbf{z}^*, \mathbf{w}^*, c/4)}} \frac{\Lambda(\mathbf{z}, \mathbf{w}) - \tilde{\Lambda}(\mathbf{z}^*, \mathbf{w}^*)}{nd} = o_P \left( \frac{\|\mathbf{z} - \mathbf{z}^*\|_{0,\sim}}{n} + \frac{\|\mathbf{w} - \mathbf{w}^*\|_{0,\sim}}{d} \right)$$

□
## B.4. Proof of Proposition 5.8 (contribution of local assignments)

**Proof.**

By Proposition 4.2, it is enough to prove that the sum is small compared to $p(\mathbf{z}^*, \mathbf{w}^*, \mathbf{x}; \boldsymbol{\theta}^*)$ on $\Omega_1$. We work conditionally on $(\mathbf{z}^*, \mathbf{w}^*) \in \mathcal{Z}_1 \times \mathcal{W}_1$. Choose $(\mathbf{z}, \mathbf{w})$ in $S(\mathbf{z}^*, \mathbf{w}^*, \tilde{c})$. This set is non empty as soon as $\min(\tilde{c}n, \tilde{c}d) > 1$.

$$ \log \left( \frac{p(\mathbf{z}, \mathbf{w}, \mathbf{x}; \boldsymbol{\theta})}{p(\mathbf{z}^*, \mathbf{w}^*, \mathbf{x}; \boldsymbol{\theta}^*)} \right) = \log \left( \frac{p(\mathbf{z}, \mathbf{w}; \boldsymbol{\theta})}{p(\mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*)} \right) + F_{nd}(\boldsymbol{\theta}, \mathbf{z}, \mathbf{w}) $$

We can assume without loss of generality that $(\mathbf{z}, \mathbf{w})$ is the representative closest to $(\mathbf{z}^*, \mathbf{w}^*)$ and denote $r_1 = \|\mathbf{z} - \mathbf{z}^*\|_0$ and $r_2 = \|\mathbf{w} - \mathbf{w}^*\|_0$. Then:

$$ \begin{aligned} F_{nd}(\boldsymbol{\theta}, \mathbf{z}, \mathbf{w}) &\le \Lambda(\mathbf{z}, \mathbf{w}) - \tilde{\Lambda}(\mathbf{z}, \mathbf{w}) + \tilde{\Lambda}(\mathbf{z}, \mathbf{w}) \\ &\le \Lambda(\mathbf{z}, \mathbf{w}) - \tilde{\Lambda}(\mathbf{z}, \mathbf{w}) - \frac{c\delta(\boldsymbol{\alpha}^*)}{4} (dr_1 + nr_2) \\ &\le -\frac{c\delta(\boldsymbol{\alpha}^*)}{4} (dr_1 + nr_2) (1 + o_P(1)) \end{aligned} $$

where the first line comes from the definition of $\Lambda$, the second line from Proposition 4.7 and the fact that $\tilde{c} < C$ and the third from Proposition 5.7 and the fact that $\tilde{c} < c/4$. Thanks to corollary D.3, we also know that:

$$ \log \left( \frac{p(\mathbf{z}, \mathbf{w}; \boldsymbol{\theta})}{p(\mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*)} \right) \le \mathcal{O}_P(1) \exp \{M_{c/4}(r_1 + r_2)\} $$

There are at most $\binom{n}{r_1}\binom{n}{r_2}g^{r_1}m^{r_2}$ assignments $(\mathbf{z}, \mathbf{w})$ at distance $r_1$ and $r_2$ of $(\mathbf{z}^*, \mathbf{w}^*)$ and each of them has at most $g^m m^m$ equivalent configurations. Therefore,

$$ \begin{aligned} & \frac{\sum_{(\mathbf{z}, \mathbf{w}) \in S(\mathbf{z}^*, \mathbf{w}^*, \tilde{c})} p(\mathbf{z}, \mathbf{w}, \mathbf{x}; \boldsymbol{\theta})}{\frac{p(\mathbf{z}^*, \mathbf{w}^*, \mathbf{x}; \boldsymbol{\theta}^*)}{(\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)}} \\ &\le \mathcal{O}_P(1) \sum_{r_1+r_2 \ge 1} \binom{n}{r_1} \binom{n}{r_2} g^{g+r_1} m^{m+r_2} \exp \left( (r_1+r_2)M_{c/4} - \frac{c\delta(\boldsymbol{\alpha}^*)}{4} (dr_1+nr_2)(1+o_P(1)) \right) \\ &= \mathcal{O}_P(1) \left( 1 + e^{\log g + M_{c/4} - d \frac{c\delta(\boldsymbol{\alpha}^*)(1+o_P(1))}{4}} \right)^n \left( 1 + e^{\log m + M_{c/4} - n \frac{c\delta(\boldsymbol{\alpha}^*)(1+o_P(1))}{4}} \right)^d - 1 \\ &\le \mathcal{O}_P(1) a_{nd} \exp(a_{nd}) \end{aligned} $$

where $a_{nd} = ne^{\log g + M_{c/4} - d \frac{c\delta(\boldsymbol{\alpha}^*)(1+o_P(1))}{4}} + de^{\log m + M_{c/4} - n \frac{c\delta(\boldsymbol{\alpha}^*)(1+o_P(1))}{4}} = o_P(1)$ as soon as $n \gg \log d$ and $d \gg \log n$.

□
## B.5. Proof of Proposition 5.9 (contribution of equivalent assignments)

### Proof.

Choose $(s, t)$ permutations of $\{1, \dots, g\}$ and $\{1, \dots, m\}$ and assume that $\mathbf{z} = \mathbf{z}^{*,s}$ and $\mathbf{w} = \mathbf{w}^{*,t}$. Then $p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \theta) = p(\mathbf{x}, \mathbf{z}^{*,s}, \mathbf{w}^{*,t}; \theta) = p(\mathbf{x}, \mathbf{z}^{*}, \mathbf{w}^{*}; \theta^{s,t})$. If furthermore $(s, t) \in \text{Sym}(\theta)$, $\theta^{s,t} = \theta$ and immediately $p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \theta) = p(\mathbf{x}, \mathbf{z}^{*}, \mathbf{w}^{*}; \theta)$. We can therefore partition the sum as

$$ \sum_{(\mathbf{z}, \mathbf{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)} p(\mathbf{x}, \mathbf{z}, \mathbf{w}; \theta) = \sum_{s,t} p(\mathbf{x}, \mathbf{z}^{*,s}, \mathbf{w}^{*,t}; \theta) \\ = \sum_{s,t} p(\mathbf{x}, \mathbf{z}^{*}, \mathbf{w}^{*}; \theta^{s,t}) \\ = \sum_{\theta' \sim \theta} \# \text{Sym}(\theta') p(\mathbf{x}, \mathbf{z}^{*}, \mathbf{w}^{*}; \theta') \\ = \# \text{Sym}(\theta) \sum_{\theta' \sim \theta} p(\mathbf{x}, \mathbf{z}^{*}, \mathbf{w}^{*}; \theta') $$

The complete likelihood $p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta)$ is a unimodal function of $\theta$ with mode located in $\hat{\theta}_c$. By consistency of $\hat{\theta}_c$, either $p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta) = o_P(p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*))$ or $p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta) = O_P(p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*))$ when $\theta$ is in a close neighborhood of $\theta^*$. In the latter case, any $\theta' \sim \theta$ other than $\theta$ is bounded away from $\theta^*$ and thus $p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta') = o_P(p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*))$. In summary,

$$ \sum_{\theta' \sim \theta} \frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta')}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*)} = \max_{\theta' \sim \theta} \frac{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta')}{p(\mathbf{x}, \mathbf{z}^*, \mathbf{w}^*; \theta^*)} (1 + o_P(1)) $$

<figure><img src="image_10.png" /></figure>

## Appendix C: Concentration for sub-exponential variables

Concentration inequalities for sub-exponential variables play a key role: in particular Proposition C.4 for global convergence and Propositions C.2 and C.3 for local convergence. We present here some properties of sub-exponential variables [14], then derives the needed concentration inequalities.

Recall first that a random variable $X$ is sub-exponential with parameters $(\tau^2, b)$ if for all $\lambda$ such that $|\lambda| \le 1/b$,

$$ \mathbb{E}[e^{\lambda(X - \mathbb{E}(X))}] \le \exp\left(\frac{\lambda^2 \tau^2}{2}\right). $$

In particular, all distributions coming from a natural exponential family are sub-exponential. Sub-exponential variables satisfy a large deviation Bernstein-type
inequality:

$$ \mathbb{P}(X - \mathbb{E}[X] \ge t) \le \begin{cases} \exp\left(-\frac{t^2}{2\tau^2}\right) & \text{if } 0 \le t \le \frac{\tau^2}{b} \\ \exp\left(-\frac{t}{2b}\right) & \text{if } t \ge \frac{\tau^2}{b} \end{cases} $$

So that

$$ \mathbb{P}(X - \mathbb{E}[X] \ge t) \le \exp\left(-\frac{t^2}{2(\tau^2 + bt)}\right) $$

### C.1. Properties

The sub-exponential property is preserved by summation and multiplication.

* If $X$ is sub-exponential with parameters $(\tau^2, b)$ and $\alpha \in \mathbb{R}$, then so is $\alpha X$ with parameters $(\alpha^2 \tau^2, ab)$

* If the $X_i$, $i = 1, \dots, n$ are sub-exponential with parameters $(\tau_i^2, b_i)$ and independent, then so is $X = X_1 + \dots + X_n$ with parameters $(\sum_i \tau_i^2, \max_i b_i)$

Moreover, Lemma C.1 defines the sub-exponential property of the absolute value of a sub-exponential variable.

**Lemma C.1.** If $X$ is a zero mean random variable, sub-exponential with parameters $(\sigma^2, b)$, then $|X|$ is sub-exponential with parameters $(8\sigma^2, 2\sqrt{2b})$.

**Proof.**

Denote $\mu = \mathbb{E}|X|$ and consider $Y = |X| - \mu$. Choose $\lambda$ such that $|\lambda| < (2\sqrt{2b})^{-1}$. We need to bound $\mathbb{E}[e^{\lambda Y}]$. Note first that $\mathbb{E}[e^{\lambda Y}] \le \mathbb{E}[e^{\lambda X}] + \mathbb{E}[e^{-\lambda X}] < +\infty$ is properly defined by sub-exponential property of $X$ and we have

$$ \mathbb{E}[e^{\lambda Y}] \le 1 + \sum_{k=2}^{+\infty} \frac{|\lambda|^k \mathbb{E}[|Y|^k]}{k!} $$

where we used the fact that $\mathbb{E}[Y] = 0$. We know bound odd moments of $|\lambda Y|$.

$$ \mathbb{E}[|\lambda Y|^{2k+1}] \le (\mathbb{E}[|\lambda Y|^{2k}]\mathbb{E}[|\lambda Y|^{2k+2}])^{1/2} \le \frac{1}{2}(\lambda^{2k}\mathbb{E}[Y^{2k}] + \lambda^{2k+2}\mathbb{E}[Y^{2k+2}]) $$

where we used first Cauchy-Schwarz and then the arithmetic-geometric mean inequality. The Taylor series expansion can thus be reduced to

$$ \begin{aligned} \mathbb{E}[e^{\lambda Y}] &\le 1 + \left(\frac{1}{2} + \frac{1}{2 \cdot 3!}\right) \mathbb{E}[Y^2]\lambda^2 + \sum_{k=2}^{+\infty} \left(\frac{1}{(2k)!} + \frac{1}{2}\left[\frac{1}{(2k-1)!} + \frac{1}{(2k+1)!}\right]\right) \lambda^{2k} \mathbb{E}[Y^{2k}] \\ &\le \sum_{k=0}^{+\infty} 2^k \frac{\lambda^{2k} \mathbb{E}[Y^{2k}]}{(2k)!} \\ &\le \sum_{k=0}^{+\infty} 2^{3k} \frac{\lambda^{2k} \mathbb{E}[X^{2k}]}{(2k)!} = \cosh(2\sqrt{2\lambda X}) = \mathbb{E}\left[\frac{e^{2\sqrt{2\lambda X}} + e^{-2\sqrt{2\lambda X}}}{2}\right] \\ &\le e^{\frac{8\lambda^2 \sigma^2}{2}} \end{aligned} $$
where we used the well-known inequality $E[|X - E[X]|^k] \le 2^k E[|X|^k]$ to substitute $2^{2k} E[X^{2k}]$ to $E[Y^{2k}]$.

<figure><img src="image_11.png" /></figure>

### C.2. Concentration inequalities

**Proposition C.2 (Maximum in $(z, w)$).** Let $(z, w)$ be a configuration and $\hat{x}_{k,l}(z, w)$ resp. $\bar{x}_{k,l}(z, w)$ be defined in Equations (3.1) and (4.1). Under the assumptions of the section 2.2, for all $\varepsilon > 0$

$$ \mathbb{P} \left( \max_{z,w} \max_{k,l} \hat{\pi}_k(z) \hat{\rho}_l(w) | \hat{x}_{k,l} - \bar{x}_{k,l}| > \varepsilon \right) \le 2g^{n+1} m^{d+1} \exp \left( - \frac{nd\varepsilon^2}{2(\bar{\sigma}^2 + \kappa^{-1}\varepsilon)} \right). \quad (C.1) $$

Additionally, the suprema over all $c/2$-regular assignments satisfies:

$$ \mathbb{P} \left( \max_{z \in \bar{Z}_1, w \in W_1} \max_{k,l} |\hat{x}_{k,l} - \bar{x}_{k,l}| > \varepsilon \right) \le 2g^{n+1} m^{d+1} \exp \left( - \frac{ndc^2\varepsilon^2}{8(\bar{\sigma}^2 + \kappa^{-1}\varepsilon)} \right). \quad (C.2) $$

Note that equations C.1 and C.2 remain valid when replacing $c/2$ by any $\tilde{c} < c/2$.

**Proof.**

The random variables $X_{ij}$ are subexponential with parameters $(\bar{\sigma}^2, 1/\kappa)$. Conditionally to $(z^*, w^*)$, $z_{+k}w_{+l}(\hat{x}_{k,l} - \bar{x}_{k,l})$ is a sum of $z_{+k}w_{+l}$ centered subexponential random variables. By Bernstein's inequality [11], we therefore have for all $t > 0$

$$ \mathbb{P}(z_{+k}w_{+l}|\hat{x}_{k,l} - \bar{x}_{k,l}| \ge t) \le 2 \exp \left( - \frac{t^2}{2(z_{+k}w_{+l}\bar{\sigma}^2 + \kappa^{-1}t)} \right) $$

In particular, if $t = ndx$,

$$ \begin{aligned} \mathbb{P}(\hat{\pi}_k(z) \hat{\rho}_l(w) | \hat{x}_{k,l} - \bar{x}_{k,l}| \ge x) &\le 2 \exp \left( - \frac{ndx^2}{2(\hat{\pi}_k(z) \hat{\rho}_l(w) \bar{\sigma}^2 + \kappa^{-1}x)} \right) \\ &\le 2 \exp \left( - \frac{ndx^2}{2(\bar{\sigma}^2 + \kappa^{-1}x)} \right) \end{aligned} $$

uniformly over $(z, w)$. Equation (C.1) then results from a union bound. Simi-
larly,

$$
\begin{align*}
\mathbb{P}(|\hat{x}_{k,\ell} - \bar{x}_{k\ell}| \ge x) &= \mathbb{P}(\hat{\pi}_k(\mathbf{z})\hat{\rho}_\ell(\mathbf{w})|\hat{x}_{k,\ell} - \bar{x}_{k\ell}| \ge \hat{\pi}_k(\mathbf{z})\hat{\rho}_\ell(\mathbf{w})x) \\
&\le 2 \exp\left(-\frac{ndx^2\hat{\pi}_k(\mathbf{z})^2\hat{\rho}_\ell(\mathbf{w})^2}{2(\hat{\pi}_k(\mathbf{z})\hat{\rho}_\ell(\mathbf{w})\bar{\sigma}^2 + \kappa^{-1}x\hat{\pi}_k(\mathbf{z})\hat{\rho}_\ell(\mathbf{w}))}\right) \\
&\le 2 \exp\left(-\frac{ndc^2x^2}{8(\bar{\sigma}^2 + \kappa^{-1}x)}\right)
\end{align*}
$$

Where the last inequality comes from the fact that c/2-regular assignments
satisfy $\hat{\pi}_k(\mathbf{z})\hat{\rho}_\ell(\mathbf{w}) \ge c^2/4$. Equation (C.2) then results from a union bound
over $\mathcal{Z}_1 \times W_1 \subset \mathcal{Z} \times W$.

<figure><img src="image_12.png" /></figure>

**Proposition C.3 (Maximum in non equivalent (z, w)).** Let $(\bar{z}, \bar{w})$ be any configuration and $(\mathbf{z}, \mathbf{w})$ the $\sim$-equivalent configuration that achieves $\|\mathbf{z} - \mathbf{z}^*\|_0 = \|\bar{z} - \mathbf{z}^*\|_0$, and $\|\mathbf{w} - \mathbf{w}^*\|_0 = \|\bar{w} - \mathbf{w}^*\|_0$, let $\hat{x}_{k\ell} = \hat{x}_{k,\ell}(\mathbf{z}, \mathbf{w})$ (resp. $\bar{x}_{k\ell}(\mathbf{z}, \mathbf{w})$) and $\hat{x}_{k\ell}^* = \hat{x}_{k,\ell}(\mathbf{z}^*, \mathbf{w}^*)$ (resp. $\bar{x}_{k\ell}^* = \bar{x}_{k,\ell}(\mathbf{z}^*, \mathbf{w}^*) = \psi'(\alpha_{k\ell}^*)$) be as defined in Equations (3.1) and (4.1). Under the assumptions $(H_1)$ to $(H_3)$, for all $\varepsilon \le \kappa\bar{\sigma}^2$,

$$
\mathbb{P} \left( \max_{(\bar{z}, \bar{w}) \sim (\mathbf{z}^*, \mathbf{w}^*)} \max_{k, \ell} \frac{nd [\hat{\pi}_k(\mathbf{z})\hat{\rho}_\ell(\mathbf{w})(\hat{x}_{k,\ell} - \bar{x}_{k\ell}) - \hat{\pi}_k(\mathbf{z}^*)\hat{\rho}_\ell(\mathbf{w}^*)(\hat{x}_{k\ell}^* - \bar{x}_{k\ell}^*)]}{n \|\mathbf{w} - \mathbf{w}^*\|_0 + d \|\mathbf{z} - \mathbf{z}^*\|_0} > \varepsilon \right) \\
\leq g^{g+1} m^{m+1} a_{nd} e^{a_{nd}}
$$

$$
\text{where } a_{nd} = nge^{-\frac{d\varepsilon^2}{\bar{\sigma}^2}} + mde^{-\frac{n\varepsilon^2}{\bar{\sigma}^2}}.
$$

**Proof.**

Denote $r_1 = \|\mathbf{z} - \mathbf{z}^*\|_0$ and $r_2 = \|\mathbf{w} - \mathbf{w}^*\|_0$. The numerator within the max in the fraction can be expanded to

$$
Z_{k\ell}(\mathbf{z}, \mathbf{w}) = \sum_{i,j} (z_{ik} w_{j\ell} - z_{ik}^* w_{j\ell}^*)(X_{ij} - \alpha_{z_{ik}^* w_{j\ell}^*})
$$

and is thus a sum of at most $N = nr_2 + dr_1$ non-null centered sub-exponential
random variables with parameters $(\bar{\sigma}^2, 1/\kappa)$. It is therefore centered sub-exponential
with parameters $(N\bar{\sigma}^2, 1/\kappa)$. By Bernstein inequality, for all $\varepsilon \le \kappa\bar{\sigma}^2$ we have

$$
\mathbb{P}(Z \ge \varepsilon(nr_2 + dr_1)) \le \exp\left(-\frac{(nr_2 + dr_1)\varepsilon^2}{2\bar{\sigma}^2}\right).
$$
There are at most $\binom{n}{r_1} g^r_1 g^g_2$ z at $||.||_{0,\sim}$ distance $r_1$ of $z^*$ and $\binom{d}{r_2} m^{r_2} m^m$ z at $||.||_{0,\sim}$ distance $r_2$ of $w^*$. An union bound shows that:

$$
\begin{align*}
& \mathbb{P}\left( \max_{(\bar{z}, \bar{w}) \sim (z^*, w^*)} \max_{k,l} \frac{Z_{kl}(\bar{z}, \bar{w})}{n \|\bar{w} - w^*\|_0 + d \|z - z^*\|_0} \ge \varepsilon \right) \\
& \le \sum_{r_1+r_2 \ge 1} \sum_{\substack{\|\bar{z} - z^*\|_{0,\sim} = r_1 \\ \|\bar{w} - w^*\|_{0,\sim} = r_2}} gm\mathbb{P}(Z_{kl}(\bar{z}, \bar{w}) \ge \varepsilon(nr_2 + dr_1)) \\
& \le \sum_{r_1+r_2 \ge 1} \binom{n}{r_1} \binom{d}{r_2} g^{g+1} m^{m+1} g^{r_1} m^{r_2} \exp\left(-(nr_2 + dr_1)\varepsilon^2/2\bar{\sigma}^2\right) \\
& = g^{g+1} m^{m+1} \sum_{r_1+r_2 \ge 1} \left(ge^{-\frac{d\varepsilon^2}{\bar{\sigma}^2}}\right)^{r_1} \left(me^{-\frac{n\varepsilon^2}{\bar{\sigma}^2}}\right)^{r_2} \le g^{g+1} m^{m+1} a_{nd} e^{a_{nd}}
\end{align*}
$$

where $a_{nd} = nge^{-\frac{d\varepsilon^2}{\bar{\sigma}^2}} + mde^{-\frac{n\varepsilon^2}{\bar{\sigma}^2}}$.

<figure><img src="image_13.png" /></figure>

**Proposition C.4 (concentration for sub-exponential).** Let $X_1, \dots, X_n$ be inde-
pendent zero mean random variables, sub-exponential with parameters $(\sigma_i^2, b_i)$.
Denote $V_0^2 = \sum_i \sigma_i^2$ and $b = \max_i b_i$. Then the random variable $Z$ defined by:

$$
Z = \sup_{\substack{\Gamma \in \mathbb{R}^n \\ \|\Gamma\|_{\infty} \le M}} \sum_i \Gamma_i X_i
$$

is also sub-exponential with parameters $(8M^2V_0^2, 2\sqrt{2Mb})$. Moreover $\mathbb{E}[Z] \leq MV_0\sqrt{n}$ so that for all $t > 0$,

$$
\mathbb{P}(Z - MV_0\sqrt{n} \ge t) \le \exp\left(-\frac{t^2}{2(8M^2V_0^2 + 2\sqrt{2Mb}t)}\right)
$$

Proof.

Note first that $Z$ can be simplified to $Z = M \sum_i |X_i|$. We just need to bound
bound $\mathbb{E}[Z]$. The rest of the proposition results from the fact that the $|X_i|$ are
subexponential $(8\sigma_i^2, 2\sqrt{2}b_i)$ by Lemma C.1 and standard properties of sums of
independent rescaled subexponential variables.

$$
\begin{align*}
\mathbb{E}[Z] &= \mathbb{E} \left[ \sup_{\|\Gamma\|_{\infty} \le M} \sum_i \Gamma_i X_i \right] = \mathbb{E} \left[ \sum_i M |X_i| \right] \le M \sum_i \sqrt{\mathbb{E}[X_i^2]} \\
&= M \sum_i \sigma_i \le M \left( \sum_i 1 \right)^{1/2} \left( \sum_i \sigma_i^2 \right)^{1/2} = MV_0 \sqrt{n}
\end{align*}
$$

using Cauchy-Schwarz.

<figure><img src="image_14.png" /></figure>
## Appendix D: Technical lemmas

Lemma D.1 is the working horse for proving Proposition 4.6. Corollary D.3 is needed for Theorem 5.8 and Lemma D.2 is an intermediate result for Corollary D.3.

### Lemma D.1.

Let $\eta$ and $\bar{\eta}$ be two matrices from $M_{g \times m}(\Theta)$ and $f: \Theta \times \Theta \to \mathbb{R}_+$ a positive function, A a (squared) confusion matrix of size g and B a (squared) confusion matrix of size m. We denote $d_{k\ell k'\ell'} = f(\eta_{k\ell}, \bar{\eta}_{k'\ell'})$. Assume that

* all the rows of $\eta$ are distinct;

* all the columns $\eta$ are distinct;

* $f(x, y) = 0 \Leftrightarrow x = y$;

* each row of A has a non zero element;

* each row of B has a non zero element;

and denote

$$ \Sigma = \sum_{kk'} \sum_{\ell\ell'} A_{kk'} B_{\ell\ell'} d_{k\ell k'\ell'} $$

Then,

$$ \Sigma = 0 \Leftrightarrow \begin{cases} A, B \text{ are permutation matrices } s, t \\ \bar{\eta} = \eta^{s,t} \text{ i.e. } \forall (k, \ell), \bar{\eta}_{k\ell} = \eta_{s(k)t(\ell)} \end{cases} $$

### Proof.

If A and B are the permutation matrices corresponding to the permutations s and t: $A_{ij} = 0$ if $i \neq s(j)$ and $B_{ij} = 0$ if $i \neq t(j)$. As each row of A contains a non zero element and as $A_{s(k)k} > 0$ (resp. $B_{s(\ell)\ell} > 0$) for all $k$ (resp. $\ell$), the following sum $\Sigma$ reduces to

$$ \Sigma = \sum_{kk'} \sum_{\ell\ell'} A_{kk'} B_{\ell\ell'} d_{k\ell k'\ell'} = \sum_{k} \sum_{\ell} A_{s(k)k} B_{t(\ell)\ell} d_{s(k)t(\ell)k\ell} $$

$\Sigma$ is null and sum of positive components, each component is null. However, all $A_{s(k)k}$ and $B_{t(\ell)\ell}$ are not null, so that for all $(k, \ell)$, $d_{s(k)t(\ell)k\ell} = 0$ and $\bar{\eta}_{k\ell} = \eta_{s(k)t(\ell)}$.

Now, if A is not a permutation matrix while $\Sigma = 0$ (the same reasoning holds for B or both). Then A owns a column $k$ that contains two non zero elements, say $A_{k_1k}$ and $A_{k_2k}$. Let $\ell \in \{1...m\}$, there exists by assumption $\ell'$ such that $B_{\ell\ell'} \neq 0$. As $\Sigma = 0$, both products $A_{k_1k}B_{\ell\ell'}d_{k_1\ell k\ell'}$ and $A_{k_2k}B_{\ell\ell'}d_{k_2\ell k\ell'}$ are zero.

$$ \left\{ \begin{array}{l} A_{k_1k} B_{\ell\ell'} d_{k_1\ell k\ell'} = 0 \\ A_{k_2k} B_{\ell\ell'} d_{k_2\ell k\ell'} = 0 \end{array} \right. \Leftrightarrow \left\{ \begin{array}{l} d_{k_1\ell k\ell'} = 0 \\ d_{k_2\ell k\ell'} = 0 \end{array} \right. \Leftrightarrow \left\{ \begin{array}{l} \eta_{k_1\ell} = \bar{\eta}_{k\ell'} \\ \eta_{k_2\ell} = \bar{\eta}_{k\ell'} \end{array} \right. \Leftrightarrow \eta_{k_1\ell} = \eta_{k_2\ell} $$

The previous equality is true for all $\ell$, thus rows $k_1$ and $k_2$ of $\eta$ are identical, and contradict the assumptions.

□
**Lemma D.2.**

Let $\mathcal{Z}_1$ be the subset of $\mathcal{Z}$ of c-regular configurations, as defined in Definition 4.1. Let $\mathbb{S}^g = \{\boldsymbol{\pi} = (\pi_1, \pi_2, \dots, \pi_g) \in [0,1]^g : \sum_{k=1}^g \pi_k = 1\}$ be the g-dimensional simplex and denote $\mathbb{S}^g_c = \mathbb{S}^g \cap [c, 1-c]^g$. Then there exists two positive constants $M_c$ and $M'_c$ such that for all $\mathbf{z}, \mathbf{z}^*$ in $\mathcal{Z}_1$ and all $\boldsymbol{\pi} \in \mathbb{S}^g_c$

$$|\log p(\mathbf{z}; \hat{\boldsymbol{\pi}}(\mathbf{z})) - \log p(\mathbf{z}^*; \hat{\boldsymbol{\pi}}(\mathbf{z}^*))| \le M_c \|\mathbf{z} - \mathbf{z}^*\|_0$$

**Proof.**

Consider the entropy map $H: \mathbb{S}^g \to \mathbb{R}$ defined as $H(\boldsymbol{\pi}) = -\sum_{k=1}^g \pi_k \log(\pi_k)$. The gradient $\nabla H$ is uniformly bounded by $\frac{M_c}{2} = \log \frac{1-c}{c}$ in $\|\cdot\|_{\infty}$-norm over $\mathbb{S}^g \cap [c, 1-c]^g$. Therefore, for all $\boldsymbol{\pi}, \boldsymbol{\pi}^* \in \mathbb{S}^g \cap [c, 1-c]^g$, we have

$$|H(\boldsymbol{\pi}) - H(\boldsymbol{\pi}^*)| \le \frac{M_c}{2} \|\boldsymbol{\pi} - \boldsymbol{\pi}^*\|_1$$

To prove the inequality, we remark that $\mathbf{z} \in \mathcal{Z}_1$ translates to $\hat{\boldsymbol{\pi}}(\mathbf{z}) \in \mathbb{S}^g \cap [c, 1-c]^g$, that $\log p(\mathbf{z}; \hat{\boldsymbol{\pi}}(\mathbf{z})) - \log p(\mathbf{z}^*; \hat{\boldsymbol{\pi}}(\mathbf{z}^*)) = n[H(\hat{\boldsymbol{\pi}}(\mathbf{z})) - H(\hat{\boldsymbol{\pi}}(\mathbf{z}^*))]$ and finally that $\|\hat{\boldsymbol{\pi}}(\mathbf{z}) - \hat{\boldsymbol{\pi}}(\mathbf{z}^*)\|_1 \le \frac{2}{n}\|\mathbf{z} - \mathbf{z}^*\|_0$.

□

**Corollary D.3.** Let $\mathbf{z}^*$ (resp. $\mathbf{w}^*$) be c/2-regular and $\mathbf{z}$ (resp. $\mathbf{w}$) at $\|\cdot\|_0$-distance c/4 of $\mathbf{z}^*$ (resp. $\mathbf{w}^*$). Then, for all $\boldsymbol{\theta} \in \boldsymbol{\Theta}$

$$\log \frac{p(\mathbf{z}, \mathbf{w}; \boldsymbol{\theta})}{p(\mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*)} \le \mathcal{O}_P(1) \exp \{M_{c/4}(\|\mathbf{z} - \mathbf{z}^*\|_0 + \|\mathbf{w} - \mathbf{w}^*\|_0)\}$$

**Proof.**

Note then that:

$$\begin{align*}
\frac{p(\mathbf{z}, \mathbf{w}; \boldsymbol{\theta})}{p(\mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\theta}^*)} &= \frac{p(\mathbf{z}, \mathbf{w}; \boldsymbol{\pi}, \boldsymbol{\rho})}{p(\mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\pi}^*, \boldsymbol{\rho}^*)} = \frac{p(\mathbf{z}, \mathbf{w}; \boldsymbol{\pi}, \boldsymbol{\rho})}{p(\mathbf{z}^*, \mathbf{w}^*; \hat{\boldsymbol{\pi}}(\mathbf{z}^*), \hat{\boldsymbol{\rho}}(\mathbf{w}^*))} \frac{p(\mathbf{z}^*, \mathbf{w}^*; \hat{\boldsymbol{\pi}}(\mathbf{z}^*), \hat{\boldsymbol{\rho}}(\mathbf{w}^*))}{p(\mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\pi}^*, \boldsymbol{\rho}^*)} \\
&\le \frac{p(\mathbf{z}, \mathbf{w}; \hat{\boldsymbol{\pi}}(\mathbf{z}), \hat{\boldsymbol{\rho}}(\mathbf{w}))}{p(\mathbf{z}^*, \mathbf{w}^*; \hat{\boldsymbol{\pi}}(\mathbf{z}^*), \hat{\boldsymbol{\rho}}(\mathbf{w}^*))} \frac{p(\mathbf{z}^*, \mathbf{w}^*; \hat{\boldsymbol{\pi}}(\mathbf{z}^*), \hat{\boldsymbol{\rho}}(\mathbf{w}^*))}{p(\mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\pi}^*, \boldsymbol{\rho}^*)} \\
&\le \exp \{M_{c/4}(\|\mathbf{z} - \mathbf{z}^*\|_0 + \|\mathbf{w} - \mathbf{w}^*\|_0)\} \times \frac{p(\mathbf{z}^*, \mathbf{w}^*; \hat{\boldsymbol{\pi}}(\mathbf{z}^*), \hat{\boldsymbol{\rho}}(\mathbf{w}^*))}{p(\mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\pi}^*, \boldsymbol{\rho}^*)} \\
&\le \mathcal{O}_P(1) \exp \{M_{c/4}(\|\mathbf{z} - \mathbf{z}^*\|_0 + \|\mathbf{w} - \mathbf{w}^*\|_0)\}
\end{align*}$$

where the first inequality comes from the definition of $\hat{\boldsymbol{\pi}}(\mathbf{z})$ and $\hat{\boldsymbol{\rho}}(\mathbf{w})$ and the second from Lemma D.2 and the fact that $\mathbf{z}^*$ and $\mathbf{z}$ (resp. $\mathbf{w}^*$ and $\mathbf{w}$) are c/4-regular. Finally, local asymptotic normality of the MLE for multinomial proportions ensures that $\frac{p(\mathbf{z}^*, \mathbf{w}^*; \hat{\boldsymbol{\pi}}(\mathbf{z}^*), \hat{\boldsymbol{\rho}}(\mathbf{w}^*))}{p(\mathbf{z}^*, \mathbf{w}^*; \boldsymbol{\pi}^*, \boldsymbol{\rho}^*)} = \mathcal{O}_P(1)$.

□
References

[1] Christophe Ambroise and Catherine Matias. New consistent and asymptotically normal parameter estimates for random-graph mixture models. *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 74(1):3–35, 2012.

[2] Peter Bickel, David Choi, Xiangyu Chang, Hai Zhang, et al. Asymptotic normality of maximum likelihood and its variational approximation for stochastic blockmodels. *The Annals of Statistics*, 41(4):1922–1943, 2013.

[3] Peter J Bickel and Aiyou Chen. A nonparametric view of network models and newman–girvan and other modularities. *Proceedings of the National Academy of Sciences*, 106(50):21068–21073, 2009.

[4] Alain Celisse, Jean-Jacques Daudin, Laurent Pierre, et al. Consistency of maximum-likelihood and variational estimators in the stochastic block model. *Electronic Journal of Statistics*, 6:1847–1899, 2012.

[5] Gérard Govaert and Mohamed Nadif. Clustering with block mixture models. *Pattern Recognition*, 36(2):463–473, 2003.

[6] Gérard Govaert and Mohamed Nadif. Block clustering with bernoulli mixture models: Comparison of different approaches. *Computational Statistics & Data Analysis*, 52(6):3233–3245, 2008.

[7] Gérard Govaert and Mohamed Nadif. Latent block model for contingency table. *Communications in Statistics—Theory and Methods*, 39(3):416–425, 2010.

[8] Gérard Govaert and Mohamed Nadif. *Co-clustering*. John Wiley & Sons, 2013.

[9] Christine Keribin, Vincent Brault, Gilles Celeux, and Gérard Govaert. Estimation and selection for the latent block model on categorical data. *Statistics and Computing*, 25(6):1201–1216, 2015.

[10] Mahendra Mariadassou and Catherine Matias. Convergence of the groups posterior distribution in latent or stochastic block models. *Bernoulli*, 21(1):537–573, 2015.

[11] Pascal Massart. *Concentration inequalities and model selection*, volume 6. Springer, 2007.

[12] J.G. Shanthikumar and U. Sumita. A central limit theorem for random sums of random variables. *Operations Research Letters*, 3(3):153 – 155, 1984.

[13] Tom A.B. Snijders and Krzysztof Nowicki. Estimation and prediction for stochastic blockmodels for graphs with latent block structure. *Journal of Classification*, 14(1):75–100, Jan 1997.

[14] Martin J Wainwright. *High-dimensional statistics: A non-asymptotic viewpoint*, volume 48. Cambridge University Press, 2019.