Network Reconstruction via the Minimum Description Length Principle

Network reconstruction is the task of determining the unseen interactions between elements of a complex system when only data on their behavior are available—usually consisting of time-series or independent samples of their states. This task is required when it is either too difficult, costly, or simply impossible to determine the individual pairwise interactions via direct measurements. This is a common situation when we want to determine the associations between species only from their co-occurrence in population samples [1], the financial market dependencies from the fluctuation of stock prices [2], protein structures from observed amino-acid contact maps [3,4], gene regulatory networks from expression data [5–7], neural connectivity from functional magnetic resonance imaging (fMRI), electroencephalography (EEG), micro electrode array (MEA), or calcium fluorescence imaging data [8–10], and epidemic contact tracing from infections [11], among many other scenarios.

In this work we frame network reconstruction as a statistical inference problem [12], where we assume that the observed data have been sampled from a generative statistical model that has a weighted network as part of its set of parameters. In this case, the reconstruction task consists of estimating these parameters from data. This generative framing contrasts with other inferential approaches based on, e.g., Granger causality [13] and transfer entropy [14] in important ways: (1) It allows for a principled framework that does not rely on ad hoc and often arbitrary confidence thresholds, (2) can reliably distinguish between direct and indirect interactions, to the extent that this information is obtainable from the data, (3) enables uncertainty quantification and comparison between alternative generative hypotheses, and (4) produces a generative model as a final result, which can be used not only for interpretation, but also to predict the outcomes of interventions and alternative conditions previously unseen in the data. These same features make this type of approach preferable also to noninferential heuristics such as correlation thresholding [15,16], known to suffer from substantial biases [12,17].

A central obstacle to the effective implementation of inferential network reconstruction is the need for statistical regularization: A network with 𝑁 nodes can in principle have up to 𝑂⁡(𝑁2) edges, and elementary maximum likelihood approaches do not allow us to distinguish between the nonexistence of an edge and the existence of an edge with a low weight magnitude, resulting in the inference of maximally dense networks, even when the true underlying network is sparse, i.e., its number of edges is far smaller than the number of possible edges. In this case, unregularized maximum likelihood would massively overfit, incorporating a substantial amount of statistical noise in its description, and thus would not meaningfully represent the true underlying system. Because of this, a robust inferential framework needs to include the ability to determine the most appropriate model complexity from data, which in this context means the right number of edges and their weight distribution, in a manner that does not overfit.

Perhaps the most common strategy to avoid overfitting is 𝐿1 regularization [18], where a penalty proportional to the sum of the weight magnitudes is subtracted from the likelihood, multiplied by a regularization parameter. This approach is a standard technique in statistics and machine learning, perhaps most well known as the central ingredient of LASSO regression [19], and has been used extensively as well for network reconstruction using multivariate Gaussian models [20–23], most famously as part of the glasso algorithm [24], but also for the inverse Ising model [25–31], among other problem instances.

The 𝐿1 regularization scheme has certain attractive properties: (1) It is easy to implement, (2) it results in the existence of exact zeros for the reconstructed weights, and (3) the 𝐿1 penalty function is convex. When combined with a likelihood function that is also convex (such as the one for the multivariate Gaussian [24] or the Ising model [29]), the latter property means that this scheme yields an overall convex objective function, thus enabling the use of very efficient algorithms available only for this simpler class of optimization problems. On the other hand, this scheme also has well-known significant limitations [18]: (1) The regularization parameter, which determines the final sparsity of the reconstruction, needs to be provided prior to inference, requiring the desired level of sparsity to be known beforehand, and (2) it causes “shrinkage” of the weights, i.e., the weights of the nonzero edges also decrease in magnitude, introducing a bias. The first disadvantage undermines this regularization approach in a central way, since it requires as an input what it should in fact provide as an output: the network sparsity. Because of this, the scheme is often used together with cross-validation [24,32,33] to determine the optimal value of the regularization parameter. But not only does this significantly increase the computational cost of the approach, but in the end also typically results in substantial overfitting, as we demonstrate in this work.

Other regularization approaches also have been proposed for network reconstruction, such as 𝛿 thresholding [34], decimation [35,36], EBIC [37], adaptive LASSO [38–40], MCP [41], SCAD [42], weight integration [43], alternative priors [44–46] with fixed number of edges [47], as well as spike-and-slab [48,49] and horseshoe [50–52] shapes. These alternatives address some limitations of 𝐿1, in particular its detrimental shrinkage properties, but they all either still require cross-validation, rely on heuristics to determine the resulting sparsity, and/or make strong assumptions on the shape of the weight distribution.

To the best of our knowledge, there is currently no regularization scheme that is simultaneously (1) principled, (2) effective at preventing overfitting, (3) algorithmically efficient, and (4) nonparametric, in particular without requiring the network sparsity to be known in advance.

In this work we fill this gap by developing a Bayesian regularization scheme that follows the minimum description length (MDL) principle [53,54] and seeks to find the optimal sparsity and weight distribution of the edges in a manner that most compresses the observed data. Instead of weight shrinkage, our approach relies on weight quantization as a means of quantifying the information necessary to encode the inferred weights to a sufficient numerical precision, and also makes no strong assumption on the shape of the weight distribution, unlike 𝐿1 and most Bayesian approaches previously proposed for this problem [43–45,47–52]. The resulting approach does not require cross-validation and is also algorithmically efficient, especially when combined with a recent algorithm that can find the best potential edge candidates to be updated in subquadratic time and in parallel [55], making it applicable for the reconstruction of networks with hundreds of thousands of nodes, or even more depending on available computing resources. Combined with the fact that it does not rely on any property of the generative model other than it being conditioned on real-valued edge weights, our approach gives us a regularization scheme that can be readily employed on a broad class of inferential reconstruction problems.

This paper is divided as follows. In Sec. II we describe our inferential framework, and in Sec. III we discuss the problems with 𝐿1 regularization. In Sec. IV we present our alternative MDL framework, and in Sec. V the inference algorithm we propose to implement it, which we then evaluate with synthetic and empirical data. In Sec. VI we conduct an empirical case study where we reconstruct the large-scale network of interactions of microbial species from empirical data and showcase how the inferred models can be used to predict the outcomes of interventions, as well as to identify “keystone” species and tipping points. We finalize in Sec. VII with a conclusion.

The general scenario for inferential network reconstruction consists of some data 𝑿 that are assumed to originate from a generative model with a likelihood

where 𝑾 ∈ℝ𝑁×𝑁 is a symmetric matrix corresponding to the weights of an undirected graph of 𝑁 nodes. In general, we are not interested in any particular a priori constraint on 𝑾, although we usually expect it to be sparse, i.e., its number of nonzero entries scales as 𝑂⁡(𝑁). In many cases, the data are represented by a 𝑁 ×𝑀 matrix of 𝑀 samples, with 𝑋𝑖⁢𝑚 being a value associated with node 𝑖 for sample 𝑚, such that

with 𝒙𝑚 being the 𝑚th column of 𝑿. Alternatively, we may have that the network generates a Markovian time series with likelihood

given some initial state 𝒙0. Other variations are possible, but for our purposes we need only to refer to a generic posterior distribution,

that we need to be able to compute up to a normalization constant. Our overall approach is designed to work independently of what kind of generative model is used, static or dynamic, with any type of likelihood, as long as it is conditioned on edge weights according to the above equation. In Appendix A we list some generative models that we will use as examples in this work, but for now it is simpler if we consider the generative model more abstractly.

Given Eq. (4) we can proceed in a variety of ways, for example, by sampling from it or computing expectations. In this work, we are interested in the maximum a posteriori (MAP) point estimate,

corresponding to the most likely weighted network that generated the data. In this setting, regularization is implemented by choosing a suitable prior 𝑃⁡(𝑾), which should, ideally, properly penalize overly complex models.

In the following we will discuss 𝐿1 regularization, and then present our alternative approach based on the MDL principle.

The most common choice to regularize Eq. (6) is independent Laplace distributions,

which result in a 𝐿1 penalty to the log-likelihood:

This choice offers some attractive properties. First, the function log 𝑃⁡(𝑾|𝜆) is convex, and therefore if log 𝑃⁡(𝑿|𝑾) is also convex with respect to 𝑾 (which is true for key problem instances such as the multivariate Gaussian or the Ising model), then Eq. (6) amounts to a convex optimization problem, which has a unique global solution with no other local maxima, and can be solved efficiently [55]. Second, this type of penalization promotes sparsity, i.e., as 𝜆 is increased, this will cause the values of 𝑊𝑖⁢𝑗 to converge in succession exactly to zero at specific values of 𝜆. This is in contrast to, e.g., 𝐿2 regularization [56], where an increase in the penalty tends to merely reduce the magnitude of the weights continuously, without making them exactly zero. In this way, the value of 𝜆 amounts to a proxy for the density of the inferred network.

However, this choice also suffers from important shortcomings. From a modeling perspective, the prior of Eq. (7) corresponds to a maximum entropy distribution conditioned on the expected mean weight magnitude, including the zeros. This amounts to a prior that is conditioned on previous knowledge about the density of the network, encoded via the parameter 𝜆. However, this information is rarely known in advance, and is almost always one of the primary objectives of inference in the first place. The likelihood of Eq. (7) is not helpful in this regard, since it is unbounded, and diverges toward log 𝑃⁡(𝑾|𝜆) →∞ for 𝑊𝑖⁢𝑗 =0 and 𝜆 →∞; therefore we cannot infer 𝜆 via MAP estimation.

The most common solution employed in practice is to select 𝜆 via 𝐾-fold cross-validation [24,32,33], i.e., by dividing the data matrix in 𝐾 disjoint subsets 𝑿 = {𝑿1,…,𝑿𝐾}, and reconstructing the network with one subset 𝑘 removed, i.e.,

and using the held-out data to compute the log-likelihood ℒ𝑘⁡(𝜆)=log 𝑃⁡(𝑿𝑘|^𝑾𝑘⁡(𝜆)). Finally, 𝜆 is chosen to maximize the mean held-out log-likelihood:

This standard technique is designed to prevent overfitting, since it relies on the predictive performance of the inferred model ^𝑾𝑘⁡(𝜆) with respect to the held-out data 𝑿𝑘, which should deteriorate as soon as it incorporates more statistical noise. Although widespread, this solution still suffers from some key limitations. First, it increases the computational cost by at least 1 or 2 orders of magnitude, since we need to perform the optimization of Eq. (9) 𝐾 times (a typical choice is 𝐾 =5) and for many values of 𝜆—usually via a combination of grid and bisection search. Furthermore, and even more importantly, this whole procedure does not sufficiently reconcile the prevention of overfitting with the promotion of sparsity, as we now demonstrate.

For our analysis it will be useful to evaluate the accuracy of a reconstruction ^𝑾 via the Jaccard similarity 𝑠⁡(𝑾,^𝑾)∈[0,1],

and likewise for the binarized adjacency matrix, 𝐴𝑖⁢𝑗⁡(𝑾)={1  if  |𝑊𝑖⁢𝑗|>0,else 0}, with 𝑠⁡(𝑨⁡(𝑾),𝑨⁡(^𝑾)).

In Fig. 1 we consider the reconstruction of Zachary’s karate club network [57], with 𝑁 =34 nodes and 𝐸 =78 edges serving as the nonzero couplings of a kinetic Ising model (see Appendix A), with each corresponding value of 𝑾 sampled from a normal distribution with mean 𝑁/2⁢𝐸 =0.22 and standard deviation 0.01. After sampling a single time series of 𝑀 =1000 transitions given a random initial state, we performed the reconstruction for a range of values of 𝜆, including also a fivefold cross-validation procedure for each value. As seen in Fig. 1(a), the increase in 𝜆 simultaneously sparsifies and shrinks the magnitude of the inferred edge weights. The maximum of the mean held-out log-likelihood obtained from cross-validation is achieved for a value 𝜆 ≈30.8, for a network with 𝐸 =278 nonzero weights, shown in Fig. 1(b). The weighted and unweighted similarities with the true network are 0.83 and 0.44, respectively, affected primarily by the abundance of spurious edges incorporated into the reconstruction. If we increase instead to 𝜆 ≈246.6, the similarities become 0.67 and 0.99, respectively, with 𝐸 =79—much closer to the true value—illustrating the inherent trade-off between weight magnitude preservation and sparsity. Although for this value of 𝜆 the binarized network is much more accurate, the average weight value is only 0.11—around half of the true value—which significantly degrades the model’s predictive performance. (We observe that without knowing the true network, it would not have been possible to determine the value of 𝜆 that yields the largest unweighted similarity, since neither the model likelihood nor the cross-validation procedure offer any indication in this regard.)

FIG. 1.

It is clear from this example that the 𝐿1 regularization obtained with the Laplace prior forces an unnecessary dichotomy between sparsity and low weight magnitudes. The detrimental bias of 𝐿1 to obtain regression coefficients is well known [42] and alternatives such as adaptive LASSO [38–40], MCP [41], and SCAD [42] exist which mitigate this problem, but they all involve ad hoc parameters that need to be determined via cross-validation, without the algorithmic attractiveness of 𝐿1, since these alternatives are no longer convex. As a result, these alternatives are substantially more cumbersome to use.

Ideally, we should be able to determine which edges are better modeled by setting them to zero, without inherently reducing the magnitude of the nonzero ones. This is an inference problem on its own, and instead of committing beforehand to which kind of weight distribution is more appropriate, a more robust approach consists in formulating an open-ended class of distributions, and learning its precise shape a posteriori from the data, following the principle of maximum parsimony. This is precisely what we achieve in the next section.

Our purpose is to implement a principled regularization scheme that does not rely on weight shrinkage, promotes sparsity, and obviates the need for cross-validation. Ideally, we would like a choice of prior 𝑃⁡(𝑾) that achieves this directly via the MAP estimation:

More precisely, we want to evoke the inherent connection between Bayesian inference and data compression by interpreting the joint likelihood 𝑃⁡(𝑿,𝑾) in terms of the description length [53,54],

where the first term in the right-hand side of the last equation corresponds to the amount of information required (e.g., in bits if base 2 is used for the logarithm) to encode the data 𝑿 when the parameters 𝑾 are known, and the second term corresponds to the amount of information required to encode the parameters 𝑾. In this setting, the prior for 𝑾 should act as a penalty counteracting the likelihood term that prevents overly complex models from being inferred. The result of the optimization of Eq. (12) would then amount to the most compressive model—i.e., with the shortest description length—striking the ideal balance between quality of fit and model complexity, extirpating as much as possible statistical noise from the model description.

However, the above equivalence with compression is not valid when the weights 𝑾 are modeled as continuous values with infinite precision. In this case, their prior 𝑃⁡(𝑾) is necessarily a probability density function, and hence −log 𝑃⁡(𝑾) cannot be interpreted as information content [58]. This is exactly the reason why we cannot determine the value of 𝜆 of the Laplace prior of Eq. (7) via MAP estimation: The most likely value is not necessarily the most compressive one. In order to properly quantify the model complexity, we need be to be explicit about the precision with which we describe the parameters 𝑾, which means that they must be discretely distributed according to a probability mass function.

With the objective of precisely quantifying the model complexity, we proceed in several steps. First, we account for the sparsity of the model by introducing an auxiliary variable 𝑨 corresponding to the binary adjacency matrix that determines the values of 𝑾 that are going to be nonzero. At first, we assume 𝑨 is sampled from a uniform distribution conditioned on the total number of edges 𝐸,

(15)

and the value of 𝐸 is itself sampled from a uniform prior 𝑃⁡(𝐸)=1/[(𝑁2)+1]—ultimately an unimportant constant added simply for completeness. This choice will allow us to benefit from sparsity, as it will lead to shorter description lengths for graphs with fewer edges. Given 𝑨, we then sample the weights 𝑾 as

where 𝑃⁡(𝑊𝑖⁢𝑗) is the probability of a nonzero weight 𝑊𝑖⁢𝑗, which is the same for all (𝑖,𝑗). Importantly, in order to account for precision, 𝑃⁡(𝑊𝑖⁢𝑗) must be a probability mass function with a support over a discrete set of values. In particular, we may, for example, consider the weights to be the outcome of a quantization procedure,

where 𝑌𝑖⁢𝑗 is some continuous auxiliary value, ⌈ ·⌋ is the round operation, and Δ represents the precision considered, such that we have a conditional probability mass 𝑃⁡(𝑊𝑖⁢𝑗|𝑌𝑖⁢𝑗,Δ) =𝛿𝑊𝑖⁢𝑗,Δ⁢⌈𝑌𝑖⁢𝑗/Δ⌋. Note that we do not assume that the true values of 𝑊𝑖⁢𝑗 are actually sampled in this way; instead this prior simply articulates the exact precision Δ with which we are performing the inference.

We can then proceed by choosing a continuous distribution for the auxiliary values 𝒀, e.g., a Laplace distribution,

in which case we obtain a probability mass,

where we have excluded the value 𝑊𝑖⁢𝑗 =0 to be consistent with our parametrization based on 𝑨. Stopping at this point we would have solved two issues: (1) Via quantization we can now properly evaluate the description length of the weights and (2) we can benefit from sparsity in a manner that is independent from the weight magnitudes. However, by choosing the Laplace prior, we remain subject to weight shrinkage, since the reduction of the description length is still conditioned on their overall magnitudes—besides the discretization, Eq. (19) still amounts to 𝐿1 regularization. Therefore, our final step is to add one more level to our Bayesian hierarchy, and try to learn the weight distribution from the data, instead of assuming that the probability decays according to the Laplace distribution. For that, we assume that the discretized weights are sampled directly from an arbitrary discrete distribution,

where 𝒛 ={𝑧𝑘} is a set of 𝐾 real values representing discrete weight categories, and 𝒑 ={𝑝𝑘} are their corresponding probabilities, with ∑𝑘𝑝𝑘 =1. From this we obtain

with 𝑚𝑘 =∑𝑖<𝑗𝐴𝑖⁢𝑗⁢𝛿𝑊𝑖⁢𝑗,𝑧𝑘 being the counts of weight category 𝑘, and with a marginal distribution

assuming a uniform prior density 𝑃⁡(𝒑) =(𝐾 −1)! over the 𝐾 simplex. The latter distribution allows for a combinatorial interpretation: It corresponds to a uniform distribution among the 𝐸 edges of the weights with fixed counts 𝒎 ={𝑚𝑘}, with probability 𝑃⁡(𝑾|𝒎)=(∏𝑘𝑚𝑘!/𝐸!)⁢∏𝑖<𝑗𝛿1−𝐴𝑖⁢𝑗𝑊𝑖⁢𝑗,0, and the uniform distribution of the counts 𝒎 themselves, with probability 𝑃⁡(𝒎|𝐸)=(𝐾+𝐸−1𝐸)−1. This interpretation allows us to fix a minor inconvenience of this model, namely that it allows for a weight category 𝑧𝑘 to end up empty with 𝑚𝑘 =0, which would require us to consider category labels that never occur. Instead, we can forbid this explicitly by using a uniform prior 𝑃⁡(𝒎|𝐸)=(𝐸−1𝐾−1)−1 over strictly nonzero counts 𝒎, leading to a slightly modified model:

Finally, we need a probability mass for the weight categories 𝒛 themselves. At this point, however, the importance of this choice is minor, since we expect 𝐾 ≪𝐸. A simple choice is to evoke the quantized Laplace once more, with

with 𝑃⁡(𝒛|𝜆,Δ,𝐾) =∏𝐾𝑘=1𝑃⁡(𝑧𝑘|𝜆,Δ). Although this choice implies that there will be a residual amount of shrinkage with respect to the weight categories, it will not dominate the regularization, unless the weight categories would otherwise diverge—something that can happen with certain kinds of graphical models such as the multivariate Gaussian when pairs of nodes have exactly the same values in 𝑿. Because of this, we found this residual shrinkage to be beneficial in such corner cases, and otherwise we do not observe it to have a noticeable effect [59].

Lastly, we need a prior over the number of weight categories 𝐾, which we assume to be uniform in the allowed range, with 𝑃⁡(𝐾|𝑨) =(1 −𝛿𝐸,0)/𝐸 +𝛿𝐸,0⁢𝛿𝐾,0. Putting it all together we have

(27)

where the remaining quantities 𝒎, 𝒛, 𝐾, 𝑨, and 𝐸 in Eq. (27) should be interpreted as being functions of 𝑾. Note that this model recovers the one using Eq. (19) when 𝐾 =𝐸, i.e., every weight belongs to its own category, up to an unimportant factor 1/(𝐸 +1)—although we expect to obtain significantly smaller description length values with Eq. (27), since it can exploit regularities in the weight distribution and determine the most appropriate precision of the weights. We emphasize that, differently from Eq. (19), in Eq. (27) the parameter Δ controls the precision of the weight categories 𝒛, not directly of the weights themselves. The precision of the weights 𝑾 is primarily represented by the actual finite categories 𝒛.

The model above is conditioned on two hyperparameters, 𝜆 and Δ, which in principle can both be optimized, assuming a sufficiently uniform prior over them (therefore, contributing to a negligible constant to the description length). In practice, the effect of such an optimization is marginal, and it is more convenient to select a small value for Δ, such as Δ =10−8, which allows us to treat the weight categories as continuous for algorithmic purposes. For 𝜆 we use an initial value of 𝜆 =1, and optimize if necessary.

We observe that we have −log 𝑃⁡(𝒛|𝜆,Δ) =𝑂⁡(𝐾⁢Δ) and −log 𝑃⁡(𝑾|𝒛,𝑨) =𝑂⁢(𝐸⁢log⁡𝐸 −𝐸⁢log⁡𝐾), and therefore we will typically have 𝐾 ≪𝐸 when the description length is minimized, in which case the prior for 𝒛 will have only a minor contribution to the overall description length. Therefore, our approach is not very sensitive to alternative choices for the weight category distribution 𝑃⁡(𝒛|𝜆,Δ). We observe also that the dominating term −log 𝑃⁡(𝑾|𝒛,𝑨) from Eq. (24) is completely invariant to bijective transformations of the weights, provided the categories 𝒛 are transformed in the same manner. Therefore, any scale or shape dependence is driven solely by the relatively minor contribution of log 𝑃⁡(𝒛|𝜆,Δ), proportional only to the number of weight categories. Since the scale parameter 𝜆 can be optimized, this means that our approach does not impose any constraints on the scale and only very weak constraints on the shape of the weights (see also Ref. [59]).

Differently from the Laplace prior of Eq. (7), Eq. (27) implements a nonparametric regularization of the weights—not because it lacks parameters, but because the shape of the distribution and number of parameters adapt to what can be statistically justified from the data. In particular, the number and position of the weight categories, as well as the number of edges, will grow or shrink depending on whether they can be used to compress the data, and provide a sufficiently parsimonious network reconstruction, without overfitting.

The prior of Eq. (27) incorporates the properties we initially desired; i.e., (1) it penalizes model complexity according to the description length, (2) it exploits and promotes sparsity, and (3) it does not rely on weight shrinkage. However, these improvements come at a cost: The prior is no longer a convex function on 𝑾. Because of this, this approach requires special algorithmic considerations, which we address in the following.

Whenever Eq. (12) corresponds to a convex optimization objective, it can be performed with relatively straighforward algorithms. The simplest one is coordinate descent [61], in which every entry in the 𝑾 matrix is optimized in sequence, while keeping the others fixed. This yields an algorithm with a complexity of 𝑂⁡(𝜏⁢𝑁2), provided the optimization of each entry can be done in time 𝑂⁡(1) (e.g., via bisection search), and 𝜏 is the number of iterations necessary for sufficient convergence. A recent significant improvement over this baseline was proposed in Ref. [55], where an algorithm was presented that solves the same optimization problem in subquadratic time, typically with a 𝑂⁡(𝑁⁢log2⁡𝑁) complexity, although the precise scaling will depend on the details of the data in general. This algorithm works as a greedy version of coordinate descent, where a subset of the entries are repeatedly proposed to be updated according to the nndescent algorithm [62]. The latter implements a stochastic search that involves keeping a list of candidate edges incident on each node in the network as an endpoint, and updating this list by investigating the list of second neighbors in the candidate network. This algorithm typically converges in log-linear time, allowing the overall reconstruction algorithm also to converge in subquadratic time, in this way converging orders of magnitude faster than the quadratic baseline of the coordinate descent algorithm.

Although these algorithms cannot be used unmodified to our problem at hand, due to its nonconvexity and the presence of additional latent variables, we can use them as building blocks. We will consider different types of updates, according to the different types of latent variable we are considering, given our initial 𝑾 (initially empty with 𝑊𝑖⁢𝑗 =0) and our set of weight categories 𝒛 (also empty initially), which are applied in arbitrary order, until the description length no longer improves.

Edge updates. We use the algorithm of Ref. [55] to find the 𝜅⁢𝑁 best zero-valued entries of 𝑾 as update candidates, usually with 𝜅 =1, according to the pairwise ranking score,

where 𝑾+ is the new matrix 𝑾 with entry 𝑊𝑖⁢𝑗 replaced by the smallest positive weight category in 𝒛, and likewise with the smallest magnitude negative value for 𝑾−. If 𝒛 is currently empty, we replace the score by the gradient magnitude:

The algorithm of Ref. [55] can obtain this set in typical time 𝑂⁡(𝑁⁢log2⁡𝑁). We call ℰ the union of this set with the currently nonzero values of 𝑾. We then visit every entry (𝑖,𝑗) in ℰ, and try to perform an update 𝑊𝑖⁢𝑗 →𝑊′𝑖⁢𝑗, accepting only if it improves the objective function, chosen uniformly at random from the following options.

Edge moves. Given the same set ℰ obtained as described above, we iterate over the nodes 𝑖 of the network with nonzero degree, and for an incident edge (𝑖,𝑗) chosen uniformly at random, and an entry (𝑖,𝑢) chosen uniformly at random from the subset of ℰ with 𝑊𝑖⁢𝑢 =0, we swap their values, i.e., 𝑊𝑖⁢𝑢 ←𝑊𝑖⁢𝑗 and 𝑊𝑖⁢𝑗 ←0, and accept if this improves the objective function.

Edge swaps. We chose uniformly at random two nonzero entries in 𝑾, (𝑖,𝑗) and (𝑢,𝑣) involving four distinct nodes, and we swap their end points, i.e., (𝑊𝑖⁢𝑣,𝑊𝑖⁢𝑗) ←(𝑊𝑖⁢𝑗,𝑊𝑖⁢𝑣) and (𝑊𝑢⁢𝑗,𝑊𝑢⁢𝑣) ←(𝑊𝑢⁢𝑣,𝑊𝑢⁢𝑗), and accept only in case of improvement.

Weight category values. We iterate over each weight category 𝑧𝑘 ∈𝒛 and perform an update 𝑧𝑘 ←𝑧′𝑘, with 𝑧′𝑘 being the result of a bisection search over the allowed values.

Weight category distribution. We consider the distribution of weight categories across the nonzero entries of 𝑾 as a clustering problem, using the merge-split algorithm described in Ref. [64], composed of the following moves, which are accepted if it improves the quality function.

The workhorse of this overall procedure is the algorithm of Ref. [55] which reduces the search for candidate updates to a subquadratic time. The remaining steps of the algorithm operate on this set with linear complexity, and therefore a single “sweep” of all latent variables can be accomplished under the same algorithmic complexity as Ref. [55], since it dominates. In practice, due to the extra latent variables and optimization steps, this algorithm finishes in time around 5–10 times longer than when using 𝐿1 regularization for a single value of 𝜆, but is typically much faster than doing cross-validation and scanning over many values of 𝜆. A reference c++ implementation of this algorithm is freely available as part of the graph-tool python library [65].

A. Assessment of the algorithm with synthetic data

In Fig. 2 we demonstrate the use of our algorithm on synthetic data generated from empirical networks. We simulate the kinetic and equilibrium Ising models [66] (see Appendix A) with true weights sampled i.i.d. from a normal distribution with mean 𝜇 =1/⟨𝑘⟩, with ⟨𝑘⟩ =2⁢𝐸/𝑁, and standard deviation 𝜎 =0.01. We compare our approach with a prior corresponding to the true distribution, i.e.,

𝑃⁡(𝑾|𝜇,𝜎,𝑨)=∏𝑖<𝑗[e−(𝑊𝑖⁢𝑗−𝜇)2/2⁢𝜎2√2⁢𝜋⁢𝜎]𝐴𝑖⁢𝑗⁢𝛿1−𝐴𝑖⁢𝑗𝑊𝑖⁢𝑗,0.

(30)

We observe that our MDL regularization performs very closely to the true prior, except for very few data, in which case the MDL tends to produce an empty network, whereas the true prior in fact overfits and yields a network with a number of edges larger than the true value. As we already explained, MAP estimation does not guarantee proper regularization when the priors correspond to probability densities—even if it happens to be the correct one. The overall tendency of the MDL approach is, by design, to err toward simpler models, instead of more complex ones, when the data are insufficient to provide a good accuracy.

FIG. 2.

As already discussed previously in Sec. III, 𝐿1 regularization with cross-validation performs substantially worse than our MDL approach, overfitting the data with networks that are orders of magnitude denser than they should be. Unlike using the true prior or MDL, for the examples of Fig. 2 the results obtained with 𝐿1 show no sign of converging asymptotically to true network as the data increase with 𝑀 →∞. This behavior further highlights the general unsuitability of the 𝐿1 approach even in situations where the data are abundant. We should also mention that the results obtained with 𝐿1 took typically 1–2 orders of magnitude longer to be computed, due to the need to optimize over the parameter 𝜆 using a bisection search.