Joint Estimation of Sparse Multilayer Networks via Graph Limits

Abstract: Network datasets in modern applications often involve multiple types of interactions occurring over a shared set of individuals. Characterizing the generating mechanisms of these interactions can be enhanced by joint modelling, as shared vertices allow layers to help explain the structure of other layers. We model multiplex observations using graph limits, called a scaled set of graphons, and develop a nonparametric joint estimator based on blockmodel approximations, termed the multi-network histogram. This nonparametric framework captures each layer's varying sparsity and connection structure, accounting for heterogeneity via shared latent variables across all layers. We establish the theoretical properties of the multi-network histogram, providing an upper bound for the weighted mean integrated squared error and deriving the optimal bandwidth that minimizes this error. By leveraging information across layers, this joint modelling achieves a reduction in error and a smaller optimal bandwidth, which enables high-resolution estimation even in sparser layers. Its usefulness is demonstrated through simulation studies and an application to socioeconomic networks in an Indian village.

• The paper introduces a multi-network histogram that jointly estimates layer-specific graphons and sparsity levels using shared latent groups and a common bandwidth.
• Theoretical results show that weighted estimation is governed mainly by denser layers, while homogeneous pooling gains an effective sample-size increase with more layers.
• Simulations and Indian village data show substantially better resolution for sparse layers, although consistency for estimated group labels and sparse-regime rates remain open problems.

Overview and motivation

This paper develops a nonparametric framework for jointly estimating the generative mechanisms of multiplex networks, in which $L$ layers of interactions are observed over an identical vertex set. The authors, Song and Olhede, model such data through a scaled set of graphons: a collection of symmetric, integrable functions $\{f^{(\ell)}\}_{\ell=1}^L$ coupled with layer-specific sparsity sequences $\{\rho_n^{(\ell)}\}$, such that $\rho_n^{(\ell)}f^{(\ell)}(x,y)$ is a valid edge probability. Latent positions $\xi_1,\dots,\xi_n$ are shared across layers, and edges are conditionally independent Bernoulli draws given these positions. This construction permits both the structure and the sparsity level to differ layer by layer — a degree of heterogeneity not accommodated by prior multilayer graphon methods, which assume either a common sparsity parameter (Chandna et al., 2020) or a common graphon (Navarro et al., 2022), and by recent multilayer neighbourhood smoothing extensions (He et al., 28 Jan 2026, Guo et al., 25 Feb 2026) that do not handle sparsity decaying at different orders.

The motivating example is the Indian village dataset of Banerjee et al., where twelve socioeconomic relations (borrowing money, lending money, giving advice, temple company, etc.) are observed over the same households. Layers vary widely in density and structure: the temple company layer is nearly empty, while visit come is the densest. Fitting a single-layer network histogram [olhede2014network] to each layer independently assigns each layer its own bandwidth and grouping, and the sparsest layers collapse into one large block, losing all structural resolution.

Methodology

The central estimator, the multi-network histogram, applies blockmodel approximation to each layer using identical group assignments and a common bandwidth $h$ across all layers. Group labels $\widehat{\mathbf{z}}$ are obtained by maximizing the sum of layer-wise profile likelihoods, and block heights are rescaled by the estimated layer sparsity $\widehat{\rho}_n^{(\ell)}$ to yield estimates of $f^{(\ell)}$. A special case, the homogeneous multi-network histogram, targets layers sharing a common $f$ with differing sparsity; it forms a weighted average of layer-wise block estimates with weights $\{f^{(\ell)}\}_{\ell=1}^L$0, derived as the minimizer of an upper bound on the integrated pointwise variance.

The authors justify the common labelling not as a convenience but as a structural requirement: under graph limits, a measure-preserving transformation must apply to every layer simultaneously, following the multiplexon framework (Ganguly et al., 8 Oct 2025). Computation proceeds via greedy label swapping initialized from spectral ordering of the densest layer, with bandwidth selected by a plug-in rule based on degree-slope estimation of the Hölder constants $\{f^{(\ell)}\}_{\ell=1}^L$1.

Theoretical results

The main results concern the oracle multi-network histogram, using group labels defined by ranks of the true latent positions, evaluated under a weighted mean integrated squared error (WMISE) with weights proportional to layer sparsities. Under Hölder-$\{f^{(\ell)}\}_{\ell=1}^L$2 smoothness of each $\{f^{(\ell)}\}_{\ell=1}^L$3 and standard bandwidth growth conditions, the paper establishes

$\{f^{(\ell)}\}_{\ell=1}^L$4

where $\{f^{(\ell)}\}_{\ell=1}^L$5 is the average layer sparsity. The optimal bandwidth is

$\{f^{(\ell)}\}_{\ell=1}^L$6

so it is governed primarily by the densest layers, through the sparsity-weighted average $\{f^{(\ell)}\}_{\ell=1}^L$7. Because the common bandwidth and shared grouping are determined jointly, sparser layers inherit a finer partition than they could support alone — this is the mechanism by which joint estimation reduces error in sparse layers. Theorem 1 of [olhede2014network] is recovered as the special case $\{f^{(\ell)}\}_{\ell=1}^L$8, $\{f^{(\ell)}\}_{\ell=1}^L$9.

For the homogeneous estimator, the MISE bound is $\{\rho_n^{(\ell)}\}$0 with optimal bandwidth $\{\rho_n^{(\ell)}\}$1, which is strictly narrower than the heterogeneous case: pooling $\{\rho_n^{(\ell)}\}$2 layers to estimate a single function increases effective sample size and permits finer resolution. The paper also analyzes an unweighted MISE criterion and shows its optimal bandwidth depends on the harmonic mean of the sparsities, i.e. it is driven by the sparsest layer; supplementary stress tests confirm that adding near-empty layers inflates the MISE-based bandwidth by 43–166% versus only 13–26% for the WMISE-based rule, with correspondingly larger and less predictable effects on dense-layer error.

Simulation evidence

Simulations over three structural scenarios (Homogeneous, Perturbation, Heterogeneous), four base graphons, three sparsity regimes, and $\{\rho_n^{(\ell)}\}$3, $\{\rho_n^{(\ell)}\}$4 compare the proposed methods against layer-wise network histogram, sort-and-smooth, USVT, and neighbourhood smoothing. Representative WMSE values (×100, mixed sparsity, $\{\rho_n^{(\ell)}\}$5, $\{\rho_n^{(\ell)}\}$6, Homogeneous scenario, function 1) illustrate the magnitude of the gains:

Method	$\{\rho_n^{(\ell)}\}$7	$\{\rho_n^{(\ell)}\}$8	$\{\rho_n^{(\ell)}\}$9
mnhist	3.537	2.888	2.412
h-mnhist	2.798	2.100	1.523
nethist	11.439	11.465	11.523
SAS	3.197	3.238	3.273
USVT	10.398	9.147	8.573
NBS	17.988	18.517	19.246

Three findings stand out. First, the proposed methods' error decreases with the number of layers, whereas competitors' errors are flat or increasing, since only the joint methods exploit cross-layer information; supplementary decomposition attributes this reduction primarily to improved label estimation rather than bandwidth selection (mnhist with estimated labels is 1.4–5× worse than with oracle labels, while the plug-in bandwidth is within a factor of 1.0–1.9 of the oracle bandwidth). Second, layer-wise results show the largest improvements in the sparsest layers — e.g. in the Homogeneous scenario with $\rho_n^{(\ell)}f^{(\ell)}(x,y)$0, sparsest-layer MSE is 12.168 for mnhist versus 64.052 for the layer-wise network histogram, and 1.523 for the homogeneous estimator. Third, the empirical decay rates match the theory in dense settings ($\rho_n^{(\ell)}f^{(\ell)}(x,y)$1) but fall short in the all-sparse regime, where observed decay is close to $\rho_n^{(\ell)}f^{(\ell)}(x,y)$2 rather than the theoretical $\rho_n^{(\ell)}f^{(\ell)}(x,y)$3 — a discrepancy the paper reports without resolving.

Application to Indian village networks

For Village 40 (231 households, 12 layers, edge densities 0.0021–0.0198), the data-driven bandwidth is $\rho_n^{(\ell)}f^{(\ell)}(x,y)$4, producing ten groups of roughly 23 households. The temple company layer — which collapses to a single block of bandwidth 231 under the single-layer method — is resolved at the joint bandwidth, revealing which household groups co-attend religious services. The estimated groups align with covariates the method never observed: groups 1 and 9 are predominantly scheduled caste, group 6 mostly general caste, and groups 3, 5, and 8 rely mainly on private electricity. Using a higher-order two-sample network test with Bonferroni correction, nine layers are found statistically indistinguishable, and the homogeneous multi-network histogram is fitted to them with bandwidth 12, yielding a finer common graphon estimate as Theorem 2 predicts.

Limitations and open questions

The paper is explicit about several constraints. The theoretical guarantees are established only for oracle group labels; no consistency result is proved for the labels estimated by the greedy algorithm. The authors note that the minimax risk for graphon estimation separates into a nonparametric term and a clustering term of order $\rho_n^{(\ell)}f^{(\ell)}(x,y)$5 in the single-layer case [gao2015rate], and that consistency of joint latent position estimation across multiple networks has only begun to be studied — in dense settings with a single common graphon (Sogan et al., 16 Mar 2026) or with an aggregated graphon via ordinal embedding — whereas the present setting with distinct graphons and diverging sparsity orders is harder and remains open. A second limitation is the two-stage layer-specific bandwidth variant, whose merging rule assumes expected degree is monotone in the latent position; this fails for graphons with non-monotone or constant degree functions, connecting to known identifiability issues in graphon models (Sogan et al., 1 Jul 2026). Third, the framework assumes identical vertex sets and no inter-layer edges; extending to perturbed or partially shared vertex sets, and to supra-adjacency formulations capturing inter-layer edges, is left as future work. Finally, the sparse-regime gap between observed ($\rho_n^{(\ell)}f^{(\ell)}(x,y)$6) and theoretical ($\rho_n^{(\ell)}f^{(\ell)}(x,y)$7) error decay is reported but unexplained.

Conclusion

This paper extends blockmodel-based nonparametric graphon estimation to multiplex networks with layer-specific sparsity and structure, via the scaled set of graphons and the multi-network histogram. Its principal contributions are a WMISE theory showing that the optimal common bandwidth is controlled by the densest layers, a homogeneous variant whose error decreases with the number of layers, and empirical demonstrations — synthetic and on Indian village data — that joint estimation materially improves resolution in sparse layers. The main theoretical gap is the absence of consistency guarantees under estimated labels, which the authors identify as the natural next step.