Analytic tractability of generalized hypergeometric network ensembles: closed-form or scalable moments and null distributions for the multivariate Wallenius adjacency distribution under arbitrary dyadic propensities

Generalized hypergeometric ensembles (gHypEG) model multi-edge networks as edge samples drawn without replacement from an urn whose colours correspond to node pairs, with odds biased by a dyadic propensity matrix; the resulting adjacency distribution is the multivariate Wallenius noncentral hypergeometric distribution. The framework enables model selection and hypothesis testing via the Mahalanobis distance between an observed network and the ensemble, but three computational and analytic gaps in the founding formulation remain unresolved: (1) the expected adjacency matrix has no closed form and is obtained only by numerically solving nonlinear defining equations, a route that comes with no accuracy, stability or complexity guarantees at the scale of large networks under arbitrary propensities; (2) the adjacency covariance matrix is available only through Monte Carlo approximation; (3) the null distribution of the Mahalanobis test statistic admits a closed-form complementary cumulative distribution only in a restrictive asymptotic regime (uniform propensities and abundant urn capacities), so p-values are otherwise obtained by extensive sampling. Solving these problems — through closed-form expressions, certified scalable algorithms, or provably accurate analytic null distributions valid in realistic finite regimes — would convert gHypEG-based network hypothesis testing from a sampling-bound enterprise into an analytically tractable one.

Background

A central task in network science is to decide whether an observed topological pattern (e.g., a community structure or a tendency for particular pairs of nodes to interact) is statistically surprising once node degrees are accounted for. This requires null models: probability distributions over networks that preserve given aggregate statistics, such as the number of edges and (exactly or in expectation) the degree sequence. The classical configuration model fixes the degree sequence exactly and is usually sampled via random rewiring; soft-constrained alternatives such as the Chung–Lu model treat edges as independent. Both routes have drawbacks: exact rewiring gives no tractable likelihood, while independent-edge models cannot keep the number of edges fixed.

The generalized hypergeometric ensemble (gHypEG), introduced by Casiraghi, Nanumyan, Scholtes and Schweitzer in 2016, interpolates between these extremes. For a directed multi-edge network (one in which a pair of nodes can be joined by more than one edge, as in interaction or mobility data) with $n$ nodes and exactly $m$ edges, one constructs an urn: for each ordered pair of nodes $(i,j)$ the urn contains $\Xi_{ij}=k_i^{\mathrm{out}}k_j^{\mathrm{in}}$ balls of a distinct colour, where $k_i^{\mathrm{out}}$ and $k_j^{\mathrm{in}}$ are the observed out- and in-degrees. Drawing $m$ balls from this urn without replacement generates the unbiased ensemble, whose adjacency distribution is the (central) multivariate hypergeometric distribution; the degrees are preserved in expectation while the number of edges is preserved exactly. Biasing the draws by a dyadic propensity matrix $\boldsymbol{\Omega}$ — each colour $(i,j)$ is drawn with relative odds $\Omega_{ij}$, so that only the ratios $\Omega_{ij}/\Omega_{kl}$ are meaningful and $\boldsymbol{\Omega}$ is defined up to a common positive factor — yields the generalized ensemble, whose adjacency matrix $\mathbf{A}$ follows the multivariate Wallenius noncentral hypergeometric distribution, i.e., the distribution of colour counts when balls are sampled without replacement with unequal odds. Arbitrary dyadic patterns (block or community structure, reciprocity effects, covariate-driven propensities) can be encoded through $\boldsymbol{\Omega}$, which makes gHypEG a flexible family of analytically specified null models with a closed-form probability mass function (up to a one-dimensional numerical integral).

Statistical testing with gHypEG uses the squared Mahalanobis distance, a multivariate generalization of the z-score measuring how many standard deviations the observed adjacency matrix $\hat{\mathbf{A}}$ lies from the ensemble mean $\langle\mathbf{A}\rangle$ relative to the covariance matrix $\boldsymbol{\Sigma}$ of the adjacency entries:

$$D_r^2(\hat{\mathbf{A}}) = (\hat{\mathbf{A}} - \langle\mathbf{A}\rangle_r)^{\mathsf{T}}\, \boldsymbol{\Sigma}_r^{-1}\,(\hat{\mathbf{A}} - \langle\mathbf{A}\rangle_r).$$

Because the number of edges is fixed — every realization satisfies $\sum_{ij}A_{ij}=m$ — the covariance matrix $\boldsymbol{\Sigma}_r$ is necessarily singular: the all-ones vector lies in its null space, so $\mathrm{rank}\,\boldsymbol{\Sigma}_r\le n^2-1$, in line with the $n^2-1$ degrees of freedom that appear below. The inverse in the quadratic form is accordingly understood as a generalized (Moore–Penrose) inverse, equivalently as the inverse of $\boldsymbol{\Sigma}_r$ restricted to the hyperplane $\{\sum_{ij}A_{ij}=m\}$; since $\boldsymbol{\Sigma}_r v=0$ implies that $v^{\mathsf{T}}\mathbf{A}$ is almost surely constant, every generalized inverse gives the same value of $D_r^2$.

A p-value for a hypothesis encoded by propensities $\boldsymbol{\Omega}_0$ is the complementary cumulative distribution function (ccdf) of $D_0^2$ under the corresponding ensemble, evaluated at the observation. Three computational obstructions arise from the routes adopted in the founding letter and persist through the framework's subsequent development; the letter itself presents them as matter-of-fact computational choices rather than as open problems. First, the mean $\langle A_{ij}\rangle$ of the multivariate Wallenius distribution is not available in closed form: it is characterized only implicitly, as the solution of a system of nonlinear equations for which the reference numerical methods are those of Fog (2008). In the definitive journal treatment that system (Corollary 5) equates the quantities $(1-\langle A_{ij}\rangle/\Xi_{ij})^{1/\Omega_{ij}}$ over all dyads and adds the constraint $\sum_{ij}\langle A_{ij}\rangle=m$, so all entries are determined by a single common constant and the whole expected adjacency matrix follows from one scalar root-find; the open issue is therefore not the number of solves, but the absence of any closed form and of published accuracy, stability and complexity guarantees for this numerical route on large networks under arbitrary $\boldsymbol{\Omega}$. Second, the covariance matrix $\boldsymbol{\Sigma}$ of the noncentral ensemble has no closed form either; the original work approximates it numerically via Fog's Monte Carlo sampling methods. Third, the ccdf of the Mahalanobis statistic has a known closed form only in a restrictive asymptotic regime — a $\chi^2$ limit with $n^2-1$ degrees of freedom that requires uniform propensities ($\Omega_{ij}=c$ for all pairs) and urn capacities much larger than the number of edges ($\Xi_{ij}\gg m$) — so in general, p-values are obtained by simulating thousands of ensemble realizations. The companion software (the ghypernet R package) still obtains the expected adjacency matrix by numerical root-finding on the defining equations and computes p-values by (parallel) Monte Carlo sampling, with a multinomial approximation used as a heuristic fallback for large networks.

Problem Statement

For a generalized hypergeometric ensemble (gHypEG) on $n$ nodes with $m$ edges — the ensemble of adjacency matrices following the multivariate Wallenius noncentral hypergeometric distribution induced by a combinatorial urn matrix $\Xi$ (with $\Xi_{ij}=k_i^{\mathrm{out}}k_j^{\mathrm{in}}$ for observed degree sequences) and an arbitrary dyadic propensity matrix $\boldsymbol{\Omega}$ (positive entries, meaningful only up to a common positive factor) — make statistical inference analytically tractable by resolving three linked gaps that the founding formulation left to numerical and Monte Carlo treatment:

1. Expected adjacency matrix. Provide closed-form expressions for the entries $\langle A_{ij}\rangle$ in terms of $\Xi$, $\boldsymbol{\Omega}$ and $m$, or alternatively an efficient, numerically stable algorithm with proven complexity and accuracy guarantees for computing $\langle\mathbf{A}\rangle$ on large networks under arbitrary $\boldsymbol{\Omega}$, replacing the current implicit characterization of the mean as the solution of the defining nonlinear system (Fog's methods; Corollary 5 of the journal treatment, in which all dyads share a single common constant), which is obtained by numerical root-finding without any published certification.

2. Covariance structure. Provide analytic formulas for the covariance matrix $\boldsymbol{\Sigma}$ of the adjacency entries (including cross-dyad covariances; $\boldsymbol{\Sigma}$ is necessarily singular, of rank at most $n^2-1$) of the noncentral ensemble, or a scalable, provably accurate numerical scheme for $\boldsymbol{\Sigma}$, replacing Monte Carlo approximation as the only available route.

3. Null distribution of the test statistic. Derive the exact distribution, or general asymptotic distributions valid beyond the regime in which a closed-form ccdf is currently known (non-uniform $\boldsymbol{\Omega}$, urn capacities not much larger than $m$), of the Mahalanobis goodness-of-fit statistic $D^2=(\mathbf{A}-\langle\mathbf{A}\rangle)^{\mathsf{T}}\boldsymbol{\Sigma}^{-1}(\mathbf{A}-\langle\mathbf{A}\rangle)$, with $\boldsymbol{\Sigma}^{-1}$ the generalized inverse described above, under a gHypEG null hypothesis — or provably accurate, computationally cheap approximations of its complementary cumulative distribution function — so that p-values for network hypothesis tests can be obtained without extensive sampling. Neither the founding letter nor its journal version specifies the complementary regime, so an answer must state the asymptotic scaling (or the finite-size conditions) in which its result holds, and the sense of its validity (limit theorem, non-asymptotic error bounds, or certified approximation).

These three components jointly serve one objective: sampling-free (or cheap and certified) statistical inference — model selection and hypothesis testing — with generalized hypergeometric ensembles at the scale of real empirical networks. An answer may resolve them through closed-form expressions, certified algorithms, or asymptotic theory rather than through unbounded Monte Carlo simulation. Throughout, 'arbitrary $\boldsymbol{\Omega}$' means any admissible propensity matrix, not only uniform or block-structured ones, and complexity or accuracy guarantees are understood per instance, with their dependence on $\boldsymbol{\Omega}$, $\Xi$, $m$ and $n$ stated explicitly.

References

1. Casiraghi, G., Nanumyan, V., Scholtes, I., Schweitzer, F. (2016). Generalized Hypergeometric Ensembles: Statistical Hypothesis Testing in Complex Networks. arXiv:1607.02441, DOI 10.48550/arXiv.1607.02441, https://arxiv.org/abs/1607.02441
2. Casiraghi, G., Nanumyan, V., Scholtes, I., Schweitzer, F. (2017). From Relational Data to Graphs: Inferring Significant Links Using Generalized Hypergeometric Ensembles. Social Informatics (SocInfo 2017), Lecture Notes in Computer Science, pp. 111-120, DOI 10.1007/978-3-319-67256-4_11, https://doi.org/10.1007/978-3-319-67256-4_11
3. Casiraghi, G., Nanumyan, V. (2021). Configuration models as an urn problem. Scientific Reports 11, 13416, DOI 10.1038/s41598-021-92519-y, https://doi.org/10.1038/s41598-021-92519-y
4. Casiraghi, G. (2021). The likelihood-ratio test for multi-edge network models. Journal of Physics: Complexity 2(3), 035012, DOI 10.1088/2632-072X/ac0493, https://doi.org/10.1088/2632-072X/ac0493
5. Casiraghi, G. (2019). The block-constrained configuration model. Applied Network Science 4, 123, DOI 10.1007/s41109-019-0241-1, https://doi.org/10.1007/s41109-019-0241-1
6. Casiraghi, G., Andres, G. (2025). Empirical networks are sparse: Enhancing multiedge models with zero-inflation. PNAS Nexus 4(1), pgaf001 (online December 2024), DOI 10.1093/pnasnexus/pgaf001, https://doi.org/10.1093/pnasnexus/pgaf001
7. Fog, A. (2008). Calculation Methods for Wallenius' Noncentral Hypergeometric Distribution. Communications in Statistics - Simulation and Computation 37(2), 258-273, DOI 10.1080/03610910701790269, https://doi.org/10.1080/03610910701790269
8. Fog, A. (2008). Sampling Methods for Wallenius' and Fisher's Noncentral Hypergeometric Distributions. Communications in Statistics - Simulation and Computation 37(2), 241-257, DOI 10.1080/03610910701790236, https://doi.org/10.1080/03610910701790236
9. Casiraghi, G., Nanumyan, V., with contributions from Brandenberger, L. and Vaccario, G. ghypernet: Fit and Simulate Generalised Hypergeometric Ensembles of Graphs. R package v1.1.2, CRAN (published 2026-02-04), DOI 10.32614/CRAN.package.ghypernet, https://cran.r-project.org/package=ghypernet
10. Haines, L. M. (2025). Who's Afraid of the Wallenius Distribution? arXiv:2511.08088, DOI 10.48550/arXiv.2511.08088, https://arxiv.org/abs/2511.08088