Renormalizable Graph Embeddings For Multi-Scale Network Reconstruction

TL;DR

This paper introduces renormalizable multi-scale graph embeddings designed for reconstructing networks at various resolutions, especially useful when the original network data is incomplete or privacy-protected, with applications in economic and trade networks.

Contribution

It proposes a novel class of scale-invariant graph embeddings that improve multi-resolution network reconstruction under data limitations, extending existing maximum-entropy methods.

Findings

Demonstrated effectiveness on economic input-output networks.

Showed improved multi-scale reconstruction accuracy.

Validated on international trade networks.

Abstract

In machine learning, graph embedding algorithms seek low-dimensional representations of the input network data, thereby allowing for downstream tasks on compressed encodings. Recently, within the framework of network renormalization, multi-scale embeddings that remain consistent under an arbitrary aggregation of nodes onto block-nodes, and consequently under an arbitrary change of resolution of the input network data, have been proposed. Here we investigate such multi-scale graph embeddings in the modified context where the input network is not entirely observable, due to data limitations or privacy constraints. This situation is typical for financial and economic networks, where connections between individual banks or firms are hidden due to confidentiality, and one has to probabilistically reconstruct the underlying network from aggregate information. We first consider state-of-the-art network reconstruction techniques based on the maximum-entropy principle, which is designed to operate optimally at a fixed resolution level. We then discuss the limitations of these methods when they are used as graph embeddings to yield predictions across different resolution levels. Finally, we propose their natural 'renormalizable' counterparts derived from the distinct principle of scale invariance, yielding consistent graph embeddings for multi-scale network reconstruction. We illustrate these methods on national economic input-output networks and on international trade networks, which can be naturally represented at multiple levels of industrial and geographic resolution, respectively.

Tables

Table 1. Table 1: Table for the ION datasets reporting, for each level ℓ \ell , the number of nodes ( N ℓ N_{\ell} ), the reduced number of nodes, the total number of links and the density.

Level $\ell$	$N_{\ell }$	$L_{\ell }$	${\rho }_{\ell }$
0	972	136887	0.29
1	647	93659	0.45
2	303	30861	0.67
3	87	3364	0.89

Equations

$${\mathcal{E}}_0:=\left\{(i_0,j_0):i_0\in {\mathcal{V}}_0,j_0\in [i_0,N_0],a_{i_0{j}_0}^{(0)}=1\right\}$$

$${\mathcal{V}}_1:=\left\{i_1:={\Omega }_0(i_0)\forall i_0\in {\mathcal{V}}_0\right\}.$$

$$a_{i_1{j}_1}^{(1)}=1-\prod _{i_0\in {\Omega }_0^{-1}(i_1),j_0\in {\Omega }_0^{-1}(j_1)}(1-a_{i_0{j}_0}^{(0)})$$

$${\mathcal{V}}_{\ell +1}:=\left\{i_{\ell +1}:={\Omega }_{\ell }(i_{\ell })\forall i_{\ell }\in {\mathcal{V}}_{\ell }\right\}.$$

$${\Omega }_{0\to \ell }:={\Omega }_{\ell }\circ \cdots \circ {\Omega }_0.$$

$$a_{i_{\ell +1}{j}_{\ell +1}}^{(\ell +1)}$$

$$=1-\prod _{i_0\in {\Omega }_{\ell \to 0}^{-1}(i_{\ell +1}),j_{\ell }\in {\Omega }_{\ell \to 0}^{-1}(j_{\ell +1})}(1-a_{i_0{j}_0}^{(0)})$$

$$a_{ij}=\{\begin{array}{ll}1& p_{ij}\\ 0& 1-p_{ij}\end{array}\forall i\le j.$$

$$p_{ij}=\frac{x_i{x}_j}{1+x_i{x}_j}$$

$$\mathcal{L}(\vec{x}|\mathbf{A})$$

$$0\stackrel{!}{=}{\partial }_i\mathcal{L}(\vec{x}|\mathbf{A})\iff ⟨k_i⟩:=\sum _{j(\ne i)}{p}_{ij}=\sum _{j(\ne i)}\frac{x_i{x}_j}{1+x_i{x}_j}\stackrel{!}{=}{k}_i,$$

$$x_i\equiv \sqrt{\delta }{\varphi }_i.$$

$$p_{ij}:=\frac{\delta {\varphi }_i{\varphi }_j}{1+\delta {\varphi }_i{\varphi }_j}$$

$$⟨L⟩:=\sum _{i<j}{p}_{ij}=\sum _{i<j}\frac{\delta {\varphi }_i{\varphi }_j}{1+\delta {\varphi }_i{\varphi }_j}\stackrel{!}{=}L.$$

$${\varphi }_i:=s_i^{undirected}:=\sum _{j(\ne i)}\frac{w_{ij}+w_{ji}}{2},$$

$$p_{ij}:=p(x_{i},x_{j},w_{i})=\begin{cases}\makebox[53.56953pt][l]{$1-e^{-x_{i},x_{j}}$}i\neq j\\ \makebox[53.56953pt][l]{$1-e^{-\frac{1}{2}x_{i}^{2}-w_{i}}$}i=j\end{cases}$$

$$⟨k_i⟩:=\sum _{j(\ne i)}{p}_{ij}\stackrel{!}{=}{k}_i.$$

$$\mathcal{L}(\vec{x}|\text{A})$$

$$=\sum _{i<j}{a}_{ij}\ln (1-e^{-x_i{x}_j})-(1-a_{ij})x_i{x}_j.$$

$${\partial }_{x_i}\mathcal{L}=\sum _{j(\ne i)}\left(\frac{a_{ij}}{p_{ij}}-1\right)x_j,$$

$$\displaystyle p_{ij}=\begin{cases}\makebox[54.95238pt][l]{$1-e^{-\delta\,\phi_{i}\,\phi_{j}}$}\qquad i\neq j\\ \makebox[54.95238pt][l]{$1-e^{-\frac{\delta}{2}\,\phi_{i}^{2}-\eta w_{i}}$}\qquad i=j\end{cases}$$

$$L\stackrel{!}{=}⟨L(\delta )⟩:=\sum _{i<j}{p}_{ij}(\delta ).$$

$$\{\begin{array}{ll}{\mathbf{P}}_{sum}^{(\ell )}& \stackrel{!}{=}{\mathbf{P}}_{cog}^{(\ell )}\\ {\mathbf{P}}^{(\ell )}& \stackrel{!}{\equiv }{\mathbf{P}}^{(m)}\end{array}$$

$$x_I=\sum _{i_0\in I}{x}_{i_0}\text{ and }{w}_I=\sum _{i_0\in I}{w}_{i_0},$$

$$p_{IJ}:=p(x_{I},x_{J},w_{I})=\begin{cases}\makebox[48.0463pt][l]{$1-e^{-x_{I}x_{J}}$}I\neq J\\ \makebox[48.0463pt][l]{$1-e^{-\frac{1}{2}x_{I}^{2}-w_{I}}$}I=J.\end{cases}$$

$$p_{IJ}=\frac{x_I{x}_J}{1+x_I{x}_J}$$

$$p_{IJ}=\frac{{\delta }^{(0)}{\varphi }_I{\varphi }_J}{1+{\delta }^{(0)}{\varphi }_I{\varphi }_J}.$$

$$⟨k_i⟩:=\sum _{j(\ne i)}{p}_{ij}\equiv f(x_i)$$

$${\sigma }_m^{(\ell )}:=\frac{⟨L_{\ell }{⟩}_m-L_{\ell }}{L_{\ell }}$$

$$-\sum _{\mathbf{G}\in \mathcal{G}}P(\mathbf{G})\ln (P(\mathbf{G}))+\sum _{i\in [1,N]}{\theta }_i\left(k_i-\sum _{\mathbf{G}\in \mathcal{G}}P(\mathbf{G})k_i(\mathbf{G})\right)$$

$$p_{ij}:=\frac{e^{-\left({\theta }_i+{\theta }_j\right)}}{1+e^{-\left({\theta }_i+{\theta }_j\right)}}$$

Full text

I Introduction

Complex networks offer a solid framework for modeling a wide range of systems that support societally relevant phenomena, from economic activities to ecological processes [1, 2, 3]. Basically, any set of pair-wise interactions can be represented via a graph, by properly defining the constituents as nodes and their connections as edges. For instance, the Input-Output network (ION) [4] represents flows (links) between sectors (nodes) of an economy, the World Trade Web (WTW) [5, 6, 7, 8] considers countries as nodes and international trade flows as edges, and food webs [9] represent flows of energy among species in an ecosystem. Network science has developed a general toolkit for analysing and modeling such diverse real-world networks, often in synergy with several field-specific disciplines. Over the last few years, these interactions have expanded to the domains of machine learning and artificial intelligence. In particular, similarities between certain classes of network models introduced in statistical physics and so-called graph embedding techniques have clearly emerged [10, 11, 12, 13]. While typically developed with different goals in mind, these techniques try to encode the observed structure of a given input network onto a low-dimensional latent space, such that e.g. the original data can be compressed, downstream operations can take place in the latent space representation, and the original system becomes endowed with a simple (random) model from which alternative but realistic instances of the network can be sampled.

Recently, an important question has been raised in the context of Network Renormalization [14], namely how to consistently transform the description of a network upon a change in the chosen resolution level – which is ultimately arbitrary or dictated by data availability, thereby not necessarily reflecting any intrinsic property of the physical system. In particular, real-space [3] (geometric) and dual-space [15] (Laplacian) approaches to network renormalization have been proposed, respectively generalizing Kadanoff’s [16] and Wilson’s [17] approaches from regular lattices to heterogeneous networks. In the framework of random graph models, network renormalization leads to the requirement of invariance of a network model with respect to the chosen resolution level [18, 14]. In particular, since the constituent nodes of a network can be aggregated arbitrarily onto super-nodes or ‘blocks’, one can formulate a rigorous requirement of invariance of the functional form of the probability distribution of a random graph model under arbitrary node aggregation [18, 19, 20, 21]. In the context of models with independent edges, this requirements leads to a unique multi-scale model (MSM) for the probability of connection between nodes and to specific rules for the renormalization of the associated model parameters [18, 19]. The consequences of this result for graph embedding techniques used in machine learning have been recently highlighted: namely, while state-of-the-art node embedding models provide no guarantee of the consistency of the vector parameters representing nodes in latent space as the network resolution is varied, the ‘multi-scale’ network model mentioned above allows for the formulation of an invariant node embedding method [13]. For the first time, this method allows for the extraction of node embedding vectors that, when a set of nodes is merged into a block, transform additively into their vector sum over the set. Convenient consequences of this result are the interpretability of the (otherwise obscure) vector sum in latent space and the possibility of connecting embeddings corresponding to different resolution levels – without the need to run the embedding algorithm anew. When applied to input-output networks defined at different taxonomic levels of the industrial classification of economic sectors, multi-scale embeddings have shown improvement over alternative approaches [13].

Building on these recent results, here we address the problem of consistent multi-scale node embeddings in a modified context, i.e. when data limitations or confidentiality imply that the input network is not entirely observable. This is the typical situation e.g. in financial and production networks, where privacy requirements make the individual connections unobservable, and only aggregate information is available. For example, accessing ‘microscopic’ firm-to-firm networks would reveal trading partners and potentially the business strategies of each actor in the economic system. As a result, publicly available network data are often coarse-grained versions of the underlying microscopic networks, where the constituent firms have been aggregated either geographically into regions or countries [22] or industrially into sectors [23]. Reconstructing the underlying firm-to-firm network from partial or aggregate data remains an open challenge, and several network reconstruction methods have been proposed for this purpose [24, 25, 26]. The importance of this challenge lies in the fact that the (reconstructed) networks are used for downstream tasks such as stress tests, measurements of systemic risk, or simulations of the propagation of shocks in the economy: if coarse-grained networks are used for these tasks, the estimated systemic risk can over- or under-estimate the risk calculated at the finer scale at which the actual stress propagation takes place [27].

Our goal is proposing node embedding methods that are both consistent upon node aggregation and operable when the input graph is not fully observable, thereby connecting graph embedding methods with network reconstruction techniques. We will initially consider state-of-the-art network reconstruction methods [24, 25, 26], investigate their (lack of) consistency under node aggregation, and finally recast these methods into a multi-scale framework [13] where the fitted model parameters renormalize consistently upon coarse-graining. In particular we start with two widely used reconstruction models: the Configuration Model (CM) [28] and the fitness-induced Configuration Model (fitnCM) [29, 30]. Both models belong to the family of Exponential Random Graphs (ERGs), whose functional form is derived from the principle of maximum entropy (i.e., minimum bias) under given constraints, which is optimal for a fixed resolution scale. Specifically, the constraint for the CM is the degree sequence (number of connections per node), while for the fitnCM, it is the total number of links. Consequently, the fitnCM is more privacy-oriented, as it depends on a single global constraint and ascribes all the node-level heterogeneity to some observed fitness values (the so-called ‘fitness ansatz’). These two models are widely applied in the field of network science [24, 25, 26], motivating our choice to use them as benchmarks. We will show that, since the CM and its fitness-based variant (fitnCM) operate at a single resolution level (just like most network models), they provide an optimal reconstruction only at that level. If nodes are grouped into blocks, the parameters of these models must be re-estimated, because the parameters obtained at different hierarchical levels ($\ell=\left\{0,1\right\}$) are not explicitly unrelated to one another. We refer to this class of models as single-scale models (SSMs). Next, adapting the multiscale approach derived in [13] for the fully observable case, we introduce the multi-scale counterparts of the above models, denoting them as the degree-corrected multi-scale model (degcMSM) and the fitnessMSM (fitnMSM) respectively. These models, which are derived from the principle of scale invariance, are first evaluated at the native scale where parameters are identified (the microscopic level) and then at coarser scales, through the renormalization (summation) rule on the parameters themselves. This evaluation allows a fair comparison between the effects of the two underlying principles (maximum Shannon Entropy and scale-invariance) from which SSMs an MSMs are derived.

It should be noted that community structures (intended as partitions of nodes into sets of densely interconnected ones) are widely encountered in various networks, including those emerging in neuroscience and social sciences [31, 32]. Consequently, the methodologies introduced in this paper are applicable to a broader range of domains than the ones (Input-Output Network and World Trade Web) analyzed here. At the same time, we stress that the concept of ‘super-nodes’ or blocks used here is more general than the one of communities (as it does not require any notion of dense subgraphs), and this extended generality is one of the main advantages of the approach we are presenting. For instance, one is often forced to change the resolution level adopted to represent the input network according to some externally imposed taxonomy of nodes (e.g. the industrial classification of firms or sectors in an economic network), irrespective of whether this taxonomy has anything in common with community structure, topological node similarity or proximity, and so on.

The rest of the paper is organized as follows. In section II, we present the multi-scale setting, namely the formal definition of a graph and the coarse-graining procedure. In section III, we introduce the Configuration Model, the fitnessCM, the degree-corrected Multi-Scale Model (MSM), the maximum-likelihood MSM and the fitnessMSM, alongside with their renormalization (summation) rule. In subsection IV.1, we describe the ING Input-Output Network and the World Trade Web datasets, as well as the coarse-graining process to generate the higher-scale representations. In section V, we present and discuss the multi-scale results, focusing on network and machine-learning performance metrics. Finally, all the technical details supporting the main text are reported in the Supplementary Material [33].

II Graph Renormalization

In this section, we formally define how to represent a graph in mathematical terms and introduce the coarse-graining procedure. For clarity, quantities associated with the microscopic level are formally denoted using the subscript $\ell=0$, while those corresponding to coarser (aggregated) levels are indexed by $\ell>0$. Those quantities may be also referred to as an $\ell$-quantity. For instance, 2-vectors and 2-nodes identify the node parameters and block-nodes at level $2$, respectively. Furthermore, when a statement holds for any level $\ell\geq 0$, we adopt the notation $\xi_{ij}:=\xi^{(\ell)}_{i_{\ell}j_{\ell}}$ to refer to pair-wise components of a generic quantity $\xi$.

Let us consider a binary undirected graph with $N_{0}$ microscopic nodes $\mathcal{V}_{0}$ (indexed as $i_{0}\in[1,N_{0}]$). Their connections (called edges or links) between nodes are defined as

[TABLE]

where $a^{(0)}_{i_{0}j_{0}}=1$ indicates the presence of an edge between nodes $i_{0}$ and $j_{0}$, whereas $a^{(0)}_{i_{0}j_{0}}=0$ otherwise. This graph can be represented by an $N_{0}\times N_{0}$ adjacency matrix $\mathbf{A}^{(0)}$ which is symmetric ($a^{(0)}_{i_{0}j_{0}}=a^{(0)}_{j_{0}i_{0}}$), due to the undirected nature of the interactions. We do not consider multiple edges between the nodes, but we do allow for self-loops, which appear along the diagonal of the adjacency matrix ($a^{(0)}_{i_{0}i_{0}}=0,1$).

II.1 Coarse-Graining

To derive the coarse-grained representation of the 0-graph, we first assume that each 0-node belongs to a single, unique block. This is achieved by defining a non-overlapping, arbitrary partition $\Omega_{0}$ that maps the 0-nodes ($\mathcal{V}_{0}$) to a new set of block-nodes

[TABLE]

Secondly, two blocks are connected if at least one edge exists between their internal nodes. This relationship formally defines the coarse-grained entries $a^{(1)}_{i_{1}j_{1}}$ as

[TABLE]

where $\Omega_{0}^{-1}(i_{1})$ denotes the set of 0-nodes mapped to the block-node $i_{1}$ (and analogously for $j_{1}$). In practice, the operation performs a logical OR operation over the potential 0-links connecting two blocks. In this context, the missing 1-links reflect the non-existing connection for all the 0-pairs, while an active 1-edge may be obtained in multiple 0-cofigurations. Furthermore, the mapping $\Omega_{0}$ is surjective but not injective, since a block-node $i_{1}$ is populated by one or more 0-nodes.

We allowed for self-interactions at every order as $i_{0}$ can be equal to $j_{0}$. The feature is consistent with the emergence of self-loops at the coarse levels, which can arise either from self-loops on constituent nodes or interactions between different internal nodes. Following these steps, we derived the coarse-grained adjacency matrix $\mathbf{A}^{(1)}$ which is binary and symmetric, alike the 0-graph, but its dimension are reduced to $N_{1}\times N_{1}$ where $N_{1}$ are the total number of blocks at level $\ell=1$.

This lumping procedure can be repeated $\ell+1(\geq 1)$ times by introducing a partition $\Omega_{\ell}$ that groups the $\ell-nodes$ into $N_{\ell+1}$ block-nodes

[TABLE]

Since these partitions $\left\{\Omega_{\ell}\right\}_{\ell\geq 0}$ are nested (non overlapping), one can compose them to obtain a direct mapping of 0-nodes to $\ell$-nodes:

[TABLE]

As a result, the $(\ell+1)$-adjacency matrix ($\mathbf{A}^{(\ell+1)}$) can be obtained from the $\mathbf{A}^{(0)}$ either by coarse-graining level-by-level (see Equation 1), or in a one-step fashion through

[TABLE]

where $\Omega^{-1}_{\ell\to 0}:=\Omega^{-1}_{0}\circ\dots\circ\Omega^{-1}_{\ell}$ is the inverse map of the direct composition Equation 2.

The nested set of partitions $\left\{\Omega_{\ell}\right\}_{\ell\geq 0}$ can be uniquely described by means by a dendrogram, as shown in [18]. A “horizontal” cut of the dendrogram yields a partition at a fixed scale. For example, in the WTW case, this corresponds to aggregate countries that are located within a specific geographical distance. On the other hand, cutting the dendrogram at different heights yields multiscale clusters, such as a state and a continent.

In many applications, the context of the problem directly suggests a partitioning of the nodes into appropriate blocks, without the need for an explicit distance matrix. A possible example is the classification of economic sectors based on the NAICS taxonomy [34]. This system assigns longer codes according to a more detailed description of the industries: for example, the two-digit code 72 corresponds to Accommodation and Food Services, and it includes Full-Service Restaurants associated with $722511$.

However, when such natural partition is available, one may rely on alternative methods to obtain a partition of the 0-nodes. These include applying any standard non-overlapping community detection algorithm, e.g. Louvain [32, 15] or even assigning 0-nodes to communities at random.

III Methodology

In the previous sections, we emphasized that a single phenomenon can be analyzed at several resolution levels. Moreover, since each element of the adjacency matrix $a_{ij}^{(\ell)}\in\left\{0,1\right\}$, every edge $(i_{\ell},j_{\ell})$ can be modelled as a Bernoulli random variable. Omitting the dependence on the level $\ell$ (see section II), the adjacency matrix pairs are modeled as

[TABLE]

The models used in this study are the fitness Configuration Model (fitnCM) [30], the Configuration Model (CM) [28] (CM), the degree-corrected MSM (degcMSM) and the maximum-likelihood MSM (maxlMSM), and the fitness multi-scale model (fitnMSM) ([18]). Among these, we selected the CM as the representative of single-scale models. To the best of our knowledge, only the MSM framework is renormalizable over arbitrary node partitions.

III.1 Configuration Model

The Configuration Model (CM) [28] estimates the probability $p_{ij}$ that two nodes $(i,j)$ are connected, assuming that the empirical degree sequence $\left\{k_{i}\right\}_{i\in[1,N]}$ is known. Unlike many other network embedding methods [35], the CM relies exclusively on node degrees to learn its parameters.

By maximizing the Lagrangian in Appendix A, one obtains the following connection probability

[TABLE]

where the parameters $x_{i}:=e^{-\theta}\geq 0$ are determined by maximizing the likelihood function, excluding self-loops:

[TABLE]

where $\vec{x}:=\left[x_{i}\right]_{i\in[1,N]}$. At the stationarity point, the parameters satisfy the following systems of equations

[TABLE]

where the initial conditions are sampled uniformly at random, as the optimization problem admits a unique solution. The optimization is performed over the structural inequivalent nodes -that is, nodes with distinct degrees (see subsubsection B.2.2). As a result, nodes sharing the same degree are assigned the same parameter $x_{i}$. This establishes a one-to-one correspondence between observed degree $k_{i}$, and parameter $x_{i}$, implying that the number of independent parameters equals the number of unique degrees in the network.

A key point to emphasize is that the parameters $\vec{x}$ are fixed such that the expected degrees (left-hand side) match the observed degrees (right-hand side). Nevertheless, the solution is equivalent to the maximum likelihood estimate Equation 7, which assumes the knowledge of the full adjacency matrix. Therefore, the degree sequence serves as sufficient statistics for the Configuration Model, enabling a privacy-oriented procedure to learn the parameters. Additional details on the derivation can be found in the appendix.

III.2 fitness Configuration Model

In many real-world context, the degree sequence is not publicly available due to privacy concerns that restrict the access to node-level information, e.g. a person or a firm, [29]. However, global network metrics, such as the total number of links, can be often acquired with fewer limitations. To adapt accordingly the CM, one assumes that each hidden parameter $x_{i}$ is proportional to a node-specific observable $\phi_{i}$ referred to as “fitness”. This assumption leads to the so-called “fitness ansatz”, namely

[TABLE]

In other words, it states that the heterogeneity of the parameters is captured by the observable $\phi_{i}$ scaled by a global parameter.

Substituting this relationship into Equation 6, the fitn Configuration Model (fitnCM) probability reads

[TABLE]

where the $\delta$ parameter is determined by requiring the expected number of links $\langle L\rangle$ to match the observed link count $L$, namely

[TABLE]

As with the CM, this constraint also corresponds to the maximum likelihood estimation of the model parameters.

The choice of the fitness observable $\phi_{i}$ depends on the specific domain. For example, we set the fitness vales as the strengths (total exchange money) of a sector, formally

[TABLE]

where the $w_{ij}$ denotes the monetary flow from node $i$ to $j$ [30]. On the other hand, for the WTW, fitness are identified with the Gross-Domestic Product of a country $i$, i.e. $\phi_{i}\stackrel{{\scriptstyle!}}{{=}}GDP_{i}$ [18, 29]. In both cases, the $\phi_{i}$ are additive as the strength, or GDP, of a group of nodes is the sum of its memebers’ strengths, or GDP.

III.3 Multi-Scale Models

A natural approach for modelling a phenomenon represented at multiple scales ($\left\{\mathbf{A}\right\}:=\left\{\mathbf{A}^{(\ell)}:\ell\geq 0\right\}$) is through renormalizable models. These models retain the same functional form throughout the levels (cfr. subsection A.1), whereas their parameters renormalize. Among others (see [3]), the ideal candidate should allow for arbitrary and heterogeneous partitions $\left\{\Omega_{\ell}\right\}_{\ell\geq 0}$ of the 0-nodes, not being restricted to a specific aggregation procedure. To address these challenges, the fitness multi-scale model (fitnMSM) was proposed in [18]. In this essay, we build upon that framework by equipping it with node parameters $\vec{x}:=\left[x_{i}\right]_{i\in[1,N]}$ (see Equation 6), which can be determined either by enforcing the degree sequence (yielding the degree-corrected MSM, or degcMSM) or maximizing the likelihood (yielding the maximum-likelihood MSM, or maxlMSM). Finally, since the CM does not account for the self-loops, we also discard them in the MSMs, we introduce in the following sections.

III.3.1 Degree-Corrected MSM

The degree-corrected MSM (degcMSM) defines the probability that two nodes are connected as

[TABLE]

where $x_{i}\in\mathbb{R}_{+},w_{i}\in\left[\left.-\frac{1}{2}x_{i}^{2},\infty\right)\right.$ are seen as the propensity of creating external connections and self-loop interactions, respectively. As for the Configuration Model, we estimated the parameters based on the degree sequence $\vec{k}:=\left[k_{1},\dots,k_{N}\right]$ (no information about the self-loops) [36]. Specifically, we forced the expected degrees to coincide with the observed one, namely

[TABLE]

Since the expected degree $\langle k_{i}\rangle$ is a one-to-one (monotonic) function of the parameter $x_{i}$, each observed degree corresponds to a unique value of $x_{i}$. That is, two different nodes sharing the same observed number of neighbors have the same parameters.

III.3.2 Maximum-Likelihood MSM

In subsection III.1, enforcing the degrees (method of moments) and maximizing the likelihood gave the same parameters $\vec{x}$ (Equation 8). However, for the MSM this is no longer true, since the MSM is not an Exponential Random Graph model. Therefore, we decided to report the expected values of the maximum-likelihood MSM (maxlMSM) for a direct comparison with the degcMSM. However, this is only a “theoretical” comparison, as one should have the observed adjacency matrix, which is not available in the contex of network reconstruction. Still, it is an intriguing result to see how the maximum-likelihood solution relate with the method of moments one.

The maxlMSM uses the same connection probability of the degcMSM (Equation 12), but the parameters are found by maximizing the likelihood

[TABLE]

where $\mathbf{A}=\left\{a_{ij}\right\}_{\left\{i\leq j\right\}}$ is the adjacency matrix of the observed network. The gradient of Equation 14 reads

[TABLE]

which vanishes at the maximum point, i.e. the optimal value for the $\vec{x}$. As described in subsubsection B.2.2, at the optimal point, nodes with identical neighbors (same neighborhood, or $N(\cdot)$) have the same parameters. For example, focusing on nodes $i,j$, $\mathcal{N}(i)\stackrel{{\scriptstyle!}}{{=}}\mathcal{N}(j)\leftrightarrow x_{i}\stackrel{{\scriptstyle!}}{{=}}x_{j}$. This implies that nodes with the same degree, which share the same parameters for the degcMSM, have different maxlMSM parameters if they are connected to different nodes. However, having the same neighbors, also implies the same parameters under the degcMSM.

Optimization Details

The initial conditions for all the models have been chosen at random from the uniform distribution $\mathcal{U}[0,1]$ and used Scipy [37] to find the maximum of the likelihood. Before setting on this, we explored ways to automatically find good initial points using Particle Swarm Optimization (PSO) [38], but the results were the same of starting with random initial conditions. Further investigations on initialization strategie are beyond the scope of this paper.

III.3.3 Fitness MSM

The fitness Multi-Scale Model (fitnMSM) [18] extends the general MSM framework (see subsection B.1) to address the issue of data confidentiality—a topic discussed in more detail in subsection III.2. As with the fitness Configuration Model (fitnCM), the fitnMSM uses the fitness ansatz (Equation 9), where a single global parameter ($\delta$) controls the total number of links $L$. Differently from the fitnCM, the fitnMSM requires the fitness to be additive, that is, the fitness of an aggregated block must equal the sum of the fitness values of its individual nodes (Equation 35). Luckily, the fitness measures previously adopted for fitnCM—such as GDP or node strength—are indeed additive. Hence, the difference between the fitnCM and fitnMSM lies in the functional form of the model, not in the nature of the fitness values.

Applying the fitness ansatz, the connection probability Equation 12 becomes

[TABLE]

where $\eta\geq 0,w_{i}\in\left.\left[-\frac{\delta}{2}\,\phi_{i}^{2},\infty\right.\right)$ control the self-loops and $\delta\geq 0$ for the total number of links. Since our goal is to compare this model with the fitnCM (not including the self-loops), the $\delta$ is estimated with the method of moments, namely

[TABLE]

As discussed for the degcMSM (subsection B.1), fitting the $\delta$ with the method of moments ($\delta_{L}$) does not necessarily return the same value as the maximum-likelihood estimate ($\delta_{\mathcal{L}}$). In both ION and WTW datasets, we observed that ( $\delta_{L}>\delta_{\mathcal{L}}$ ), meaning that the likelihood-based model underestimates the actual number of links. Therefore, we used $\delta_{L}$ (method of moments), as it fixes the empirical link count as the fitnCM.

III.4 Renormalized Models

In this section, we introduce the renormalization procedure applied to the previously discussed models. The multi-scale models (MSMs) are inherently renormalizable by design (see subsection III.3). In contrast, for the CM and fitnCM, renormalizability is not inherent and must be imposed externally (see subsection A.1).

III.4.1 Summed MSM

This class of models is defined by enforcing the following two conditions, namely

[TABLE]

$\forall\ell\geq 0,m\geq 0$, which are jointly referred to as the scale-invariance (or self-consistency) properties. They require that the renormalized probability $\mathbf{P}_{sum}^{(\ell)}$ equals the coarse-grained $\mathbf{P}_{cog}^{(\ell)}$ with the same functional form ([18], subsection B.1) for all the scales.

We characterize the fitted vectors at finest level [math] as $\left\{x_{i_{0}},w_{i_{0}}\right\}_{i\in[1,N_{0}]}$. Moreover, we define $\left\{x_{I},w_{I}\right\}_{I\in[1,N_{\ell}]}$ as the block-parameters at the level $\ell>0$ where the block $I$ corresponds to a group of original nodes via $I:=\Omega_{0\to\ell}(i_{0})$. Assuming that the block-node parameters are obtained by summing the lower level ones, namely

[TABLE]

the Equation 19 is automatically satisfied due to the model’s functional form Equation 12 (see subsection B.1). In other words, this additive coarse-graining rule ensures that the structure of the model remains unchanged across all levels.

As a result, the connection probability between block nodes ( $I$) and ( $J$ ) at level ( $\ell$ ) can be computed by simply inserting the summed parameters from Equation 20 into the original MSM probability function:

[TABLE]

III.4.2 Summed CM

Once the microscopic $\vec{x}^{(0)}:=\left\{x_{i_{0}}\right\}_{i_{0}\in[1,N_{0}]}$ have been estimated, the Configuration Model (CM) lacks a built-in renormalization procedure. As such, two possible strategies can be followed: “coarse-graining” the microscopic probabilities (Equation 34) or by refitting the parameters at that level $\ell$. Both solutions suffer from a similar drawback due to the lack of scale-invariance of the CM (cfr. Equation 19). Specifically, the first method yields a coarse-grained probability matrix $\mathbf{P}_{\text{cog}}^{(\ell)}$ that can no longer be expressed in the original CM form using group-level parameters. Hence, the resulting model ceases to be a proper Configuration Model (see subsection A.1). On the other hand, refitting the CM at each level $\ell$ leads to $\ell$-vectors that are unrelated with the $\ell-1$ parameters. From the model’s perspective, each level arises from a distinct generative process, thereby losing connection of the parameters across scales.

To address these drawbacks, we used the same renormalization (summation) scheme as in Equation 20 which leads to the connection probability among (I,J) blocks:

[TABLE]

This procedure satisfies the second condition in Equation 19 by preserving the functional form of the CM at the block level. While the first condition (equivalence between renormalized and coarse-grained probabilities) is still not functionally enforced, we proceed pragmatically to determine at which extent this affects higher-level expectations.

III.4.3 Fixed (renormalized) fitnCM

Similarly to the CM, the fitnCM lacks scale-invariance properties. Therefore, we implement a renormalization scheme whereby $\delta^{(\ell)}\stackrel{{\scriptstyle!}}{{=}}\delta^{(0)}$ (the latter parameter was fixed at the microscopic level [math]), while the block-fitness $I$ is obtained by summing the inner fitness $\left\{\phi_{i}\right\}_{\left\{i\in I\right\}}$ (see Appendix A). Note that this approach is equivalent of applying the sum rule of the CM. Finally, the higher-level connection probability among (I,J) blocks is the following:

[TABLE]

IV Data

IV.1 ING Input-Output Network

ING Bank N.V. regularly reports the economic transactions of all ING clients for different years. We focused on firm-to-firm payments in the year 2022 between ING bank accounts, filtering out the transactions involving individual or non-Dutch clients, flows to or from non-ING account and the internal transfers within the same firm, i.e. self-payments. Since ING is the largest bank in The Netherlands [39], the dataset covers a major portion of the national market. More precisely, we selected the year 2022 to simplify the numerical calculations and to avoid distortions in the data from the aftermath of the COVID-19 pandemic. However, the procedure may be easily extended for other time intervals, such as multi-year spans or quarterly snapshots.

At the firm-to-firm (f2f) resolution, the dataset consists of $N_{f2f}\approx 3.4\cdot 10^{5}\text{ nodes},L_{f2f}\approx 4\cdot 10^{6}\text{ links}$, yielding a density $\rho_{f2f}\approx 3.5\cdot 10^{-5}$. This classifies the network as large and sparse. To analyze the data at a higher level, we aggregated the firms according to their NAICS (North American Industry Classification System) codes, and set an edge among two sectors if there was at least one payment among firms belonging to one of them (see Equation 3). Then, to focus on the production structure of this graph, we excluded sectors such as “Public Administration” (92), “Finance and Insurance” (52), “Management of Companies and Enterprises”/“Holdings” (55) as they are not directly related with a product/service. In particular, “Public Administration” fluxes include taxes and fees; “Finance and Insurance”, loans, that are not part of the supply chain of a product; and “Management of Companies and Enterprises”/“Holdings” include business entities that own shares in multiple firm. Lastly, we mapped, for simplicity, the 6-digits NAICS codes to integers, i.e. $111110\to 1,111120\to 2,\cdots$.

Our work is the first application of the MSM to the multiscale structure derived from the ION. Therefore, we focused on the economic relationships among the sectors, discarding both the directionality and the volume of monetary flow. In this context, each edge is undirected (reciprocated) and binary. For example111The $i,j$ refers to the observed [math]-nodes, namely $i:=i_{0},j:=j_{0}$, if $w_{ij}$ refers to the total amount of money sent from $i\textnormal{ to }j$ at level [math], we reciprocated the transaction by setting $w_{ij}\to w^{{}^{\prime}}_{ij}:=\frac{w_{ij}+w_{ji}}{2}$ [5]. Furthermore, an edge was regarded as active ($a_{ij}=a_{ji}=1$) if $w^{{}^{\prime}}_{ij}>0$ or absent otherwise. This results in a symmetrical and binary edge, whenever there was at least one flow between two sectors. At this stage, the ION consists of $N_{0}=972$ sectors and $L_{0}=1.4\cdot 10^{5}$ links, yielding a network density of $\rho_{0}\approx 0.29$, roughly 4 order of magnitudes higher than $\rho_{f2f}$.

IV.2 Coarse-Graining The ION

In the previous section, we constructed the binary undirected adjacency matrix $\mathbf{A}^{(\ell=0)}$ representing the interactions among the 6-digits sectors (0-nodes). Here, we describe the coarse-graining procedure to obtain the multi-scale unfolding of $\mathbf{A}^{(0)}$.

At first, we grouped together the 0-nodes that share the same first $6-\ell$ digits, with $\ell\in[0,4]$. For instance, the sectors $111191,111199$ are merged in the same block $11119$ starting from $\ell=1$.

Secondly, we set an edge among two $\ell$-nodes if there existed at least one connection between any of their constituent 0-nodes (see Equation 4). Our procedure uses a one-step approach to easily generate every level directly from the $\mathbf{A}^{(0)}$, without passing through the intermediate scales. The latter procedure is also allowed by the MSM.

Note that from $\ell=4$ onward, the coarse-grained graphs become fully-connected, limiting a meaningful statistical modelling up to $\ell=3$. For technical details on the evolution of the number of nodes $N_{\ell},$ links $L_{\ell}$, and density $\rho_{\ell}:=\frac{2L_{\ell}}{N_{\ell}\left(N_{\ell}-1\right)}$, refer to Table 2.

IV.3 World Trade Web

As a second application, we analyzed the World Trade Web (WTW) tracked by the Gleditsch dataset [22], which reports bilateral trade flows (imports and exports) for every country. We focused on the year 2000 (the most recent one), and excluded the missing states in the BACI-CEPII GeoDist [40], whose distances are required for the coarse-graining procedure. After this filtering step, the graph has $N_{0}=185$ 0-nodes. Although we illustrate the method for the year 2000, the methodology applies to any other period of time.

For every pair $(i,j)$, the dataset reports the exported volume $w_{ij}$ and the imported $\omega_{ji}$ one (all in USA dollars). Since the reported export from $i$ to $j$ ($w_{ij}$) generally differ from the import from $j$ to $i$ ($\omega_{ij}$), we redefined $w_{ij}\leftarrow\frac{w_{ij}+\omega_{ji}}{2}$ [41] as the averaged amount of trade from $i\textnormal{ to }j$.

Successively, as for the ION, we symmetrized the transaction matrix by mapping every weight $w_{ij}$ to the average flow between the two directions, namely $w_{ij}\to w^{{}^{\prime}}_{ij}:=\frac{1}{2}\left(w_{ij}+w_{ji}\right)$, i.e. the average flow between the two directions. By renaming $w^{{}^{\prime}}_{ij}$ as $w_{ij}$, the resulting weights are symmetric, i.e. $w_{ij}=w_{ji}$. Given the high link reciprocity of the transaction reported in the Gleditsch dataset [42], this undirected representation is a sound approximation. Then, we binarized the import-export matrix to construct the adjacency matrix $\mathbf{A}:=\Theta(\textbf{W})$ where the $\Theta(\cdot)$ is the Heaviside function.

IV.4 Coarse-Graining of WTW

To coarse-grain the level [math], we used the inter-country geographical distances [40] to iteratively merge closer nodes via single-linkage agglomerative clustering as in [18]. This procedure returns a dendrogram whose leaves are the 0-nodes, whose branching points are the block-countries, and whose branch heights equal the single-linkage distances between sub-clusters. By cutting this dendrogram at 18 successive levels $\ell\in[0,17]$ we obtain partitions $\Omega_{\ell},\ell\geq 0$, such that the number of block-countries are approximately $N_{\ell}:=N_{0}-10\cdot\ell$. For $\ell\geq 7$ the coarse-grained graph become fully connected. Therefore, a meaningful statistical modelling is possible up to $\ell=6$. See Table 2 for the evolution of the number of nodes, links, and network density.

V Results and Discussions

In this section, we apply the previously introduced models to the ION datasets. The corresponding results for the WTW dataset are provided in the Appendix.

V.1 Scale-Invariance Evidence

Figure 2 illustrates the property of scale-invariance by plotting the left-hand side (LHS) of the equation in Equation 36 against its right-hand side (RHS). By construction, both the fitnMSM and degcMSM models perfectly recover the identity line, while the maxlMSM is not depicted since it uses the same probability as the degcMSM. In contrast, the summed versions of the fitnCM and CM models consistently underestimate their coarse-grained counterparts, forming a characteristic knife-like distribution that is typical of their functional form. This finding differs from that of [13], which reported a different relationship for the LogisticPCA model [12]. A simple example demonstrating why the fitnCM and CM models lack scale-invariance is provided in Appendix A.

V.2 Summed against Fitted parameters

In the remainder of this paper, we evaluate how well the observed network measurements $Y_{i_{\ell}}$ align with the corresponding model predictions at different aggregation levels $\ell$. Specifically, we compare them to the expected values obtained using the summed probabilities $\langle Y_{i_{\ell}}\rangle_{\text{sum}}$ and the re-fitted probabilities $\langle Y_{i_{\ell}}\rangle_{\text{fit}}$.

In Figure 3, we illustrate the relationship between the summed parameters $x_{I}:=\sum_{i_{0}\in I}x_{i_{0}},I:=i_{2}$, and the fitted parameters $\hat{x}_{I}$ at aggregation level 2 for three models: the CM, degcMSM, and maxlMSM. The dashed identity line represents perfect agreement between the two.

For the CM (leftmost panel), the plot shows the characteristic “cosine-like” curve, indicating that low values of $x_{I}$ (from summation) tend to overestimate the fitted $\hat{x}_{I}$, while high values underestimate it. This bias stems from imposing a summation-based renormalization, despite the fact that the CM model is not scale-invariant and lacks a natural rule for renormalization (see Appendix A). The panels on the right correspond to the MSMs (degcMSM and maxlMSM), where the summed and fitted parameters shows a better agreement at level 2 (see also 3b and 3c), indicating better consistency under the renormalization scheme.

This discrepancy can be attributed to the differences in functional forms between the models, as shown in Figure 11. In particular, the MSM connection probability satisfies $p^{\text{MSM}}(x)>p^{\text{CM}}(x);\forall x$. As a consequence, to reproduce the same low (or high) observed degrees, the CM typically requires smaller (or larger) fitted values $\hat{x}$ compared to the MSM. Additionally, since the CM function increases more gradually, i.e. it saturates more slowly, the optimizer is more sensitive to parameter variation, which allows it to reach lower values. As a result, the CM parameter range spans nearly two orders of magnitude more than the MSM’s.

Interestingly, the close alignment between the summed and fitted parameters in the MSMs was not guaranteed a priori, although it was expected due to the multi-scale design. For a more in-depth discussion on this topic, refer to the following subsection subsubsection V.2.1.

V.2.1 Interlude: Scale-Invariance Implications

The principle of scale-invariance, as defined in Equation 33, requires that the network ensembles generated from summed and coarse-grained probabilities are statistically equivalent. This does not guarantee that the predictions derived from a summed model will match those from a re-fitted one. Therefore, true deviation from scale-invariance should be only assessed using the Equation 36 shown in Figure 2, since other metrics can lead to misleading results that depend on the specific scale or quantity.

Nevertheless, we perform this broader comparisons to investigate the practical consequences of the scale-invariance principle, which is a crucial for real-world applications. For deeper insights on how the partitions affect this devation, refer to the companion paper [13].

V.3 Fitted VS Exogenous Variables

The CM, degcMSM, and maxlMSM determine their parameters based on the network topology. For example, the CM constrains the degree sequence on average. Inspired by previous studies [29, 30], here we examine the relationship between the topologically-inferred parameters and the exogenous variables $f$, referred to as fitnesses, that characterize the underlying system222In the context of the WTW, GDP has previously been shown to play a dominant role [29]. A similar fitness-based approach was later applied to the inter-firm network (IFN) in [30]..

In Figure 4, we illustrated this relationship by plotting the node-specific parameters $x_{i}$ from the fitted and summed models on the y-axis, whereas the node strengths on the x-axis, measured at levels $\ell=0$ and $\ell=2$. To ensure additivity, at any resolution level, the strength of node $i$ is defined as: $s_{i}:=\sum_{j}w_{ij}$ which includes self-loops (i.e., monetary flows from a node to itself). Remember that we excluded the self-loops at the firm level, prior to coarse-graining, to prevent that self-payments (potentially arising from internal accounting transactions) could distort the modeling of the production network. To standardize the measurements, both parameters and fitnesses are normalized by their respective averages. Hence, the identity line in the figure represents a perfect agreement between fitness and model parameters up to a scaling factor. By construction, fitness models lie exactly on this line (with only the fitnMSM shown in the plot for clarity). Each model is represented using a distinct color, while crosses and circles indicate fitted and summed parameters, respectively, at the given resolution level.

V.3.1 Level 0

At level 0 (see 4a), node strengths closely approximate the parameters of the Configuration Model (CM), while they tend to underestimate those of the MSMs, especially for nodes with low values. As a result, the fitnCM is expected to yield better local-level predictions compared to the MSMs (refer to 7a).

As described in section III, CM and degcMSM assign the same fitted parameter to all nodes with equal degree. This behavior becomes particularly noticeable for the lowest parameter values (related with the minimum degree), where the one parameter value corresponds to various strengths.

In the case of the maxlMSM, nodes with the same neighbors receive the same parameter. The plateau observed for low-regime weights, however, arises from very small variations in the fitted parameters (on the order of ( $\sim 10^{-3}$)), rather than from identical connectivity patterns, such as all nodes linking to the same hubs.

All models, except the maxlMSM333For the ION, the signed relative error of the predicted number of links is very low, approximately $\frac{L_{\text{maxlMSM}}-L}{L}\approx 10^{-3}$. This holds true also in the context of the WTW., exactly match the observed total number of links at the fitting level. However, there is no straightforward relation between the node strengths and the expected network statistics, since strengths can either overestimate or underestimate the parameters.

For instance, the expected degree is a non-linear monotonic function $f$ of the parameter $x_{i}$, namely

[TABLE]

where $f$ depends on all the parameters but $x_{i}$. Thus, even if $s_{i}<x_{i}^{\text{degcMSM}}$, it does not necessarily imply that $\langle k_{i}\rangle{\text{fitnMSM}}<\langle k_{i}\rangle_{\text{degcMSM}}$.

For a more detailed analysis of node-level quantities, we refer the reader to the following sections.

V.3.2 Level 2

4b illustrates the behavior of the summed variables $x_{I}$ with respect to the fitted values $\hat{x}_{I}$ for the CM, degcMSM, and maxlMSM. As expected, the discrepancy in the CM is larger than that observed in the two MSMs. This is because the CM is a single-scale model that does not support renormalization by summing the internal weights, or by any other procedure, as discussed in subsection A.1.

In contrast, the parameters of the MSMs vary within a much narrower range, roughly within one order of magnitude. Interestingly, however, these models are not explicitly designed to ensure that the summed values $x_{I}$ match the fitted values $\hat{x}_{I}$. Instead, they are constructed to remain self-consistent as discussed in detail in subsubsection V.2.1. Additionally, when degree-corrected models are applied, for a single value of the fitted weight $\hat{x}_{I}$, there can be several distinct values of the corresponding summed parameter. This reflects the variability across internal strength distributions, of nodes with the same degree.

Another key observation is that most of the strength values tend to overestimate the MSM parameters, and this bias persists even when moving to coarser levels of aggregation. Surprisingly, the summed parameters in the CM, despite its lack of scale-invariance, are not drastically different from the fitted ones—highlighting an unexpected level of alignment.

V.4 Expected Number of Links

In Figure 5, we report the signed relative error $\sigma_{m}$ (on the y-axis) for each model $m$, across different resolution levels $\ell$ (on the x-axis). The signed relative error (SRE) is defined as

[TABLE]

where $\langle L_{\ell}\rangle_{m}$ is the expected number of links by the model $m$ and $L_{\ell}$ the corresponding empirical value at level $\ell$. The reconstruction models (all but the maxlMSM) were fitted to retrieve at least the total number of links at the finest level. As a consequence, $\sigma_{m}^{(0)}$ vanishes for all models except the maxlMSM, which achieves a near-zero error $\sigma_{\text{maxlMSM}}^{(0)}\approx 10^{-6}$. For higher levels of aggregation, the expected number of links is determined using the models’ summed probabilities, as defined in Equation 21 and discussed in subsection A.2.

When averaging SRE across all levels, the fitnMSM results as the most accurate model overall, followed in order by degcMSM, maxlMSM, fitnCM, and finally the CM. However, if we focus on the absolute relative error, the CM and degcMSM yield similar performance up to the final level.

More broadly, while some metrics may suggest that the CM exhibits scale-invariance, this behavior only holds up to a certain degree of resolution, beyond which its single-scale nature becomes evident. Consequently, testing for scale-invariance through a specific metric should be done with care, as discussed in subsubsection V.2.1.

V.5 Network Measurements

In Figure 6, we show how three node metrics—Degree (DEG), Average Nearest Neighbor Degree (ANND), and Clustering Coefficient (CC)—evolve across different levels of resolution. Each panel illustrates the predictions from the empirical network $\mathbf{A}^{(\ell)}$, the summed model $\mathbf{P}{\text{sum}}^{(\ell)}$, and the refitted model $\mathbf{P}{\text{fit}}^{(\ell)}$. Small inset plots display the agreement between the summed and refitted probabilities, with point-colors representing the log-density $\rho$: from blue (low density, $\rho=0$) to red (high density, $\rho=1$).

The first set of plots (6a) illustrates the behavior of the Configuration Model (CM), while the second (6b) for the degree-constrained MSM (degcMSM). Similar plots were produced for fitnCM, fitnMSM, and maxlMSM, but they did not provide depper insights than the ones appreciable from these two representative cases.

Focusing on the ANND and CC, the CM consistenty underestimates the observed metrics as the aggregation levels increase. In contrast, the degcMSM predictions remain better aligned with the empirical measurements at all scales. Interestingly, the degcMSM seems to better approximate the observed metrics also at the fitting level [math] (see also 7a) as the CM tends to overestimate ANND and CC for high-degree nodes.

These results are also reflected in the inset plots. At level 2, for example, the CM’s summed probabilities $p^{\text{CM}}{\text{sum}}$ tend to cluster near 0 and 1, underestimating the more widely distributed refitted values $p^{\text{CM}}{\text{fit}}$. On the other hand, MSM summed probabilities $p^{\text{MSM}}_{\text{sum}}$ are closely aligned with their refitted counterparts, scattered around the identity line. This means that at coarser levels, CM predictions tend to underestimate network statistics even more than MSM models do.

V.6 Reconstruction Accuracy

7a shows the Reconstruction Accuracies (RA) for the ION dataset [42]. In simple terms, RA represents the probability that an observed node-level metric $Y_{i_{\ell}}$ falls within the uncertainty range predicted by the model. Put differently, it measures how well a model captures a given metric $Y$ at resolution level $\ell$. We denote this value as $RA^{(\ell)}_{m}\left(Y\right)$ and formally define it in Equation 65.

Since all models are fitted at the microscopic level $\ell=0$, most of the highest RA values are also found at this scale. However, some exceptions occur, where the best performance appears at higher aggregation levels. A notable example is the fitnCM model, which achieves its highest RA for the clustering coefficient (CC) at level 3. This level is especially peculiar, as it is the final one before the network becomes fully connected. Consequently, the structure is much more heterogeneous than in previous levels. Additionally, the number of nodes at this stage is the smallest among all layers, meaning that the accurate prediction of even a single node contributes more significantly to the overall score.

Although in 6a we showed that the Configuration Model (CM) consistently underestimates metrics like DEG, ANND, and CC, the corresponding RAs do not drop to zero. This is because the model’s predicted uncertainty range is wide enough to include the observed values. In contrast, the degcMSM and maxlMSM—which are scale-consistent by design—achieve high RA values at all levels. The only significant exception occurs for the average nearest neighbor degree (ANND) at level 2. Here, the empirical ANND is not well captured by the MSMs. However, this appears to be a second-order effect, as we still observe strong agreement between the summed and fitted parameters (see Figure 3) and high RA values for the degrees at the same level. As discussed in subsubsection V.2.1, this is not problematic: the principle of scale-invariance in MSMs enforces self-consistency across levels, rather than exact matching of observed values at each individual resolution (see also Equation 36).

Among all models, the degree-constrained MSM (degcMSM) delivers the best performance at the fitted level $\ell=0$, outperforming the CM as well. This outcome is somewhat surprising, as previous studies showed that logistic PCA—a Shannon Entropy-based, single-scale model—gave better results than MSMs at this resolution [12]. Although the objectives of that work (network modeling) differ from ours (network reconstruction), CM was still expected to be the top performer at the microscopic scale. Yet, our results do not confirm that assumption.

Another interesting result is the performance of the maximum likelihood MSM (maxlMSM). Its reconstruction accuracies are nearly the same as those of degcMSM, suggesting that constraining degree sequences and maximizing likelihood can yield comparable results—a pattern also observed in Figure 3.

7a also includes, for completeness, the reconstruction performance of the fitness-based models. Although neither fitnCM nor fitnMSM performs well in capturing local metrics, fitnCM outperforms fitnMSM. This result aligns with Figure 4, which shows a better agreement between node strengths and topological parameters in the fitnCM case.

That said, when focusing on global metrics—such as the total number of links $L$—the fitnMSM consistently captures $L$ across scales (see Figure 5), even though we fitted the unique parameter at level [math]. This behavior aligns with their objective, i.e. capturing global metrics, where local fluctuations are average out, rather than node-level details.

V.7 ROC-PR Curves

In 7b, we reported the Area Under the Curve of the Receiver Operating Characteristic (ROC) and the Precision-Recall (PR) curves, defined in Appendix C. The CM exhibits values that are often indistinguishable from its local multi-scale variants, making it hard to find the single-scale or scale-invariant one. Moreover, the fitnCM consistently outperforms the fitnMSM, as it benefits from a more effective fitness ansatz (see Figure 4 and 7a). From these metrics, it becomes even more difficult to identify which models are scale-invariant.

Interestingly, all AUC-PR scores tend to (approximately) linearly increase with the resolution scale (see 7d). This occurs even though both true positives and predicted positives decrease with scale; their ratio, however, increases. This trend was somewhat expected: as the network becomes denser at coarser scales, the probability of correctly predicting a link among the expected ones increases—leading to a higher precision score, i.e., a larger value of $PPV$.

A similar trend appears also in AUC-ROC scores. However, at intermediate level $1$, we observe a noticeable drop in AUC-ROC scores, as the models exhibit decreased sensitivity to true positives in absolute terms. In other words, while the models improve in precision—correctly identifying positives among those they predict—they become less effective when performance is measured against the set of all true positives (recall).

The key takeaway is that evaluation highlights the relative strength of local models compared to global ones. Novertheless, it remains difficult to differentiate single-scale from multi-scale models, as their AUC performances are often coinciding. Most importantly, across both evaluation metrics and all resolution levels, fitnCM consistently outperforms fitnMSM that we know it shouldn’t be the case based on the results from Figure 5.

V.8 Expected Number of Triangles

In Figure 8, we report the observed and expected triangle densities (TriDens), defined as $\triangle(k_{i}\geq c)/N$ where the subgraph $\mathbf{G}{{i:k_{i}\geq c}}$ is obtained by extracting the nodes whose observed degree ${k_{i}}{i\in[1,N]}$ is greater than or equal to a threshold ( c ). Here, ( N ) refers to the number of nodes at the given resolution level [12].

The panel 8a shows results at level $\ell=0$, where all models were originally fitted, whereas 8b presents results at level $\ell=2$ for the summed models.

By construction, as $c$ increases, the difference among the observed and expected densities vanishes until it coincides as $c=N-1$, i.e. for the whole graph. At very high thresholds $c\lessapprox N-1$, only the highest-degree nodes (hubs) remain. However, these nodes tend to be weakly connected to each other, as the average nearest neighbor degree (ANND) decreases with increasing degree (see Figure 6). As a result, no triangle is closed among the “richest” nodes.

This outcome is not in contradiction with prior findings such as [10, 11], since we are analyzing the expected triangle densities rather than computing their exact values.

Our results show that, even with a scalar parameter (lowest-rank embedding), local models can faithfully reproduce triangle densities. In contrast, global fitness-based models have trouble replicating triangle patterns, as they are designed to capture aggregate properties of the network—such as $\sum_{i}\triangle(k_{i}\geq c)/N$—rather than node-level structural details. For this reason, results from the fitness models are included only for the sake of completeness.

VI Conclusions

The strength of graph-based models lies in their flexibility to represent different types of interactions by appropriately defining nodes and edges. By aggregating low-level (microscopic) nodes into larger block-nodes, a coarse-grained network is obtained—a simplified version of the original graph. Repeating this procedure across multiple levels produces a hierarchical, multi-scale unfolding of the observed system. We applied this approach to both the ION and WTW datasets (see subsection IV.1), showing how a single underlying generative process can be represented consistently across different resolutions.

Most existing network reconstruction techniques [26] are designed to operate at a single level of resolution, estimating link probabilities based solely on the data at that scale. As a result, they often overlook the importance of model consistency across scales.

To highlight this limitation, we compared a widely used single-scale model—the Configuration Model (CM)—with our proposed multi-scale model (see section III). For completeness, we also evaluated the performance of the fitness-based Configuration Model (fitnCM), the fitness-based MSM (fitnMSM), and the maximum-likelihood estimation of MSM parameters (maxlMSM). The key difference between these models lie in scale-invariance: MSMs are inherently renormalizable, aligning with the (hidden) generative process, whereas the CM are not. Specifically, the block-node parameters are defined as the sum of the parameters of their constituent nodes (see Equation 20), allowing for a principled interpretation of weights aggregation across scales. In contrast, CM-based models lack this scale consistency of its parameters and must be refitted independently at each level. More in general, from the single-scale model perspecive, each aggregated level is treated as an entirely new and unrelated structure—as if generated by a different underlying process—even though, in our case, all levels derive from the same multi-resolution aggregation of the original network.

As detailed in section V, the CM tends to perform well only at the resolution at which it is fitted, but struggles to generalize when applied to coarser levels. In contrast, our multi-scale model provides a consistent structural representation across resolutions, requiring only knowledge of how microscopic nodes are grouped into larger blocks. This advantage is most evident at higher aggregation levels, such as the final layer of the ION.

Together, these findings emphasize the importance of adopting a multi-scale perspective. This approach provides a robust and flexible framework for understanding complex systems by allowing analysis at various resolutions, depending on the granularity of the phenomenon under investigation.

VII DATA AND CODE AVAILABILITY

The dataset on sector transactions (ION) used in the paper is highly confidential and cannot be disclosed. However, the world trade is a long-standing public dataset, one may find at Gleditsch. Lastly, the code is freely available as the multi-scale-node-embeddings python package.

VIII ACKNOWLEDGMENTS

We thank ING Bank N.V. for their support and active collaboration. A special thanks to the whole DataScience team at ING Bank for their advice that helped shape this research. This work is supported by the projects “Reconstruction, Resilience and Recovery of Socio-Economic Networks” (RECON-NET EP FAIR 005 - PE0000013 “FAIR” - PNRR M4C2 Investment 1.3) and “SoBigData.it - Strengthening the Italian RI for Social Mining and Big Data Analytics” (Grant IR0000013 n. 3264, 28/12/2021 -https://pnrr.sobigdata.it/), financed by the European Union - NextGenerationEU. This publication is part of the project “Network renormalization: from theoretical physics to the resilience of societies” with file number NWA.1418.24.029 of the research programme NWA L3 - Innovative projects within routes 2024, which is (partly) financed by the Dutch Research Council (NWO) under the grant https://doi.org/10.61686/AOIJP05368.

Appendix A Configuration Model

The configuration model is a generative random graph model which is “maximally random” out of the node-degrees which are enforced on average. Its functional form is obtained by finding the stationarity point of the Lagrangian

[TABLE]

with respect to the probability distribution $P$. The notation used above defines $\mathcal{G}$ as the ensemble graphs G sampled from the distribution $P(\textbf{G})$ and $\theta_{i}$ is the i-th Lagrange parameter connected to the observed degree $k_{i}$. Moreover, the first term of the summation provides the “maximal uncertainty”, namely the maximum of Shannon Entropy, whereas the second imposes the expected degrees (on average) to be equal to the observed ones.

The stationarity (functional) point for Equation 26 reads

[TABLE]

where the Lagrangian parameters $\theta\in(-\infty,\infty)$ are fixed via maximum likelihood - the least unbiased functional for these kinds of models [43]. Due to exponential form, this class of models is generally addressed in Physics as Exponential Random Graph models (ERG). In our context, the Equation 6 may describe either the connection probability of a link or the likelihood that a link-fermion particle has energy $x_{i}x_{j}$ [28]. Indeed, since the links are either present or not (0 or 1 in the adjacency matrix), they must be described by the same statistics of a fermion particle.

Referring to Equation 7, the hessian of the likelihood is positive defined [36] and, therefore, the stationarity point (Equation 8) is the unique maximum of the likelihood. Moreover, it requires only the knowledge of the observed degrees -the observed quantity that one assumed to know from the beginning. Hence, this method is the “most unbiased” one [43] -no more information needed- where the degrees stand as the sufficient statistics $\left\{k_{i}\right\}_{i\in[1,N]}$ to enforce the $\vec{x}$. In general, fixing the first moment of the probability distribution -method of moments- requires less information with respect to the knowledge of the full adjacency matrix -the input of the maximum likelihood principle. However, this is enough for the CM.

A.1 Inconsistency of the CM: a trivial example

By referring to Figure 9, the microscopic network of 3 nodes $0,1,2$ merges into the coarse-grained network with blocks $A,B$, containing respectively the nodes $0,1$ and the node $2$, namely $A=\Omega(0)=\Omega(1),B=\Omega(2)$.

Assuming having fixed the CM parameters at the bottom level, the upper-level model is obtained by means of the RHS of the Equation 36. For example, the probability among $A,B$ block-nodes reads

[TABLE]

Hence, the functional dependence of $P_{cog}$ is different from the actual CM one Equation 6 since there is no renormalization rule that can aggregate the parameters in such a way to recover the Equation 6. In turns, this implies that the CM is not scale-invariant. Indeed, the LHS of the Equation 36 is

[TABLE]

which by no means equal to $P_{cog}$. Therefore, the scale-invariant principle has to be enforced as a grounding premises, changing the model functional form accordingly. Therefore, the model is not renormalizable.

A.2 Summed (renormalized) CM

To enforce the self-consistency (scale-invariance), i.e. the activation function to be the same at every scale, we impose the block-parameters are the sum of their members parameters $x_{A}:=(x_{0}+x_{1})\textnormal{ and }x_{B}:=x_{2}$ as for the MSM. With this choice, the summed probability is

[TABLE]

where $x_{I}:=\sum_{i\in I}x_{i},x_{J}:=\sum_{j\in J}x_{j}$. Therefore, the comparison between the CM and the MSM would be about how well their expected measurements can predict the observed ones.

Appendix B Multi-Scale Model

B.1 Model Derivation

By referring to the full derivation of the multi-scale model (MSM) [18], we report here only the crucial passages for our treatment. The observed network is assumed to lay at the lowest level [math], whereas the higher scales $\ell\geq 1$ identify the coarser representations obtained by aggregating the level [math] graph. However, the treatment would hold for every pair $m(\geq 0),\ell\geq m$ as pointed out in [18], to which we refer for a more rigorous derivation.

Let us consider the binary undirected adjacency matrix $\mathbf{A}^{(\ell=0)}$ describing the microscopic interactions among the 0-nodes. To obtain the coarser graphs, one needs a hierarchical and non-overlapping partition of the microscopic nodes $\left\{\mathbf{\Omega}_{\ell}\right\}_{\ell\geq 0}$ into block nodes. More concretely, the block-nodes $I:=i_{\ell}$, at level $\ell+1$, is defined as $I:=\mathbf{\Omega}_{0\to\ell-1}(i_{0})$ reached after either $\ell$ subsequent partitions or, equivalently, after one bigger aggregation $\Omega_{0\to\ell-1}:=\Omega_{\ell-1}\circ\dots\circ\Omega_{0}$. Once the groups $I$ have been identified, one assigns a link among the blocks via

[TABLE]

where $a_{IJ}:=a^{(\ell)}_{i_{\ell+1}j_{\ell+1}},a_{i_{0}j_{0}}:=a^{(0)}_{i_{0}j_{0}}$ and $J:=j_{\ell}=\mathbf{\Omega}_{0\to\ell-1}(j_{0})$. In other words, a link among the block-nodes exists if there is at least one micro-link among their members. In turns, by iterating the procedure $\ell$ times, it is possible to create $\mathbf{A}^{(\ell)}$, the low-resolution graph describing the original phenomenon the $\ell$th scale.

In order to model this architecture, one needs several assumptions. The first one requires that the MSM must describe the microscopic matrix $\textbf{A}^{(0)}$, namely \textbf{A}^{(0)}\sim P^{(0)}\big{(}\mathbf{A}^{(0)},\mathcal{X}^{(0)}\big{)} subject to P^{(0)}\big{(}\cdot,\mathcal{X}^{(0)}\big{)}\stackrel{{\scriptstyle!}}{{=}}(P^{(0)})^{T}\big{(}\cdot,\mathcal{X}^{(0)}\big{)} and \sum_{\textbf{S}\in\mathscr{A}^{(0)}}P^{(0)}\big{(}\textbf{S},\mathcal{X}^{(0)}\big{)}\stackrel{{\scriptstyle!}}{{=}}1 where $\mathscr{A}^{(0)}$ is the ensemble of all the binary symmetric graphs with $N_{0}$ nodes.

In line with [18], we further assume that the parameters

[TABLE]

are given by the product between $\left\{x_{i_{0}}\right\}_{i_{0}\in[1,N_{0}]}$ with additional node-wise parameters $\left\{w_{i_{0}}\right\}_{i_{0}\in[1,N_{0}]}$ only active in the self-loop part ($i_{0}=j_{0}$). In particular, $x_{i_{0}}$ encodes the capability of node $i_{0}$ to connect to other nodes; whereas $w_{i_{0}}$ its propensity for a self-interaction. Since the principles leading to an edge are different to the self-loop ones, we explicitly introduced the $w$ parameters. The $\frac{1}{2}$ in the self-loop part is introduced to avoid the double-counting of the self-interaction. We discarded a possible dependence of $\mathcal{X}^{(0)}$ on a dyadic $d_{i_{0}j_{0}}$ factor and the higher order terms. For further details we refer to [18].

In the following, we will use the notation P^{(0)}\big{(}\mathbf{A}^{(0)},\mathcal{X}^{(0)}\big{)}:=P^{(0)}\big{(}\mathbf{A}^{(0)},\mathbf{X}_{w}^{(0)}\big{)} where $\mathbf{X}^{(0)}_{w}:=[x^{(0)},\vec{w}^{(0)}],x^{(0)}:=\left[x_{1_{0}},\dots,x_{N_{0}}\right]^{T}\in N_{0}\times\mathbb{R}_{\geq 0},\vec{w}^{(0)}:=\left\{w_{i_{0}}\right\}_{i_{0}\in[1,N_{0}]}\in N_{0}\times[-\frac{1}{2}\left(\vec{x}_{i_{0}}\right)^{2},\infty)$ are the fitted parameters at level [math].

At the lower resolution [math], the generated ensemble $\mathscr{A}^{(0)}$ contains multiple configurations that, after coarse-graining, lead to the observed macroscopic $\mathbf{A}^{(\ell)}$, i.e. $\{\mathbf{A}^{(0)}\}\xrightarrow{\mathbf{\Omega_{0\to\ell-1}}}\mathbf{A}^{(\ell)}$ (for a nice visualization of this technical point see [18]). In turns, this induces the probability of observing $\mathbf{A}^{(\ell)}$ as

[TABLE]

The scale-invariance principle requires that the functional form of the MSM has to be independent of the chosen scale, i.e. P_{\ell}\big{(}\cdot,\cdot)\stackrel{{\scriptstyle!}}{{=}}P_{0}\big{(}\cdot,\cdot)\quad\forall\ell\geq 0, and that the model can generate the $\ell$-graphs in two equivalent ways, either hierarchically or directly. The hierarchical procedure (see Equation 31) prescribes to generate the [math]-graph ensemble with probability P\big{(}\mathbf{A}^{(0)},\mathbf{X}^{(0)}_{w}\big{)} and, then, coarse-graining them $\ell$ times via the partitions $\{\mathbf{\Omega}_{k}\}^{\ell-1}_{k=0}$. The direct method, instead, requires to calculate the renormalize parameters $\widetilde{\mathbf{X}}_{w}^{(\ell)}$ and, then, directly model $\mathbf{A}^{(\ell)}$ via P\big{(}\mathbf{A}^{(\ell)},\widetilde{\mathbf{X}}_{w}^{(\ell)}). Imposing both requirements

[TABLE]

where we have defined the LHS and RHS of the first equation as $\mathbf{P}_{sum}^{(\ell)}\textnormal{ and }\mathbf{P}_{cog}^{(\ell)}$ respectively. In other words, the form of the $P(\cdot,\cdot)$ will depend on the scale $\ell$ only through the renormalized parameters $\widetilde{\mathbf{X}}_{w}^{(\ell)}$. By assuming that the links are statistically independent, the previous equation yields

[TABLE]

where $p_{IJ}:=p_{IJ}(\widetilde{\mathbf{X}}_{w}^{(\ell)})$ and

[TABLE]

depend, respectively, on the renormalized and the fitted parameters at level $\ell=0$. The interpretation of Equation 33: the probability $p_{IJ}$ that two blocks $I,J$ are connected is given by the probability that there is at least one link among their microscopic members $i_{0}\in I,j_{0}\in J$ (see Equation 29).

The missing “sum” superscript over $p_{IJ}$ in Equation 33 reflects the scale-invariant nature of functional form, whereas the subscript $IJ$ refers to the $x_{I},x_{J}$ parameters that are, indeed, varying at coarser resolutions under renormalization.

By taking the logarithm of both sides of the Equation 33, the only functional form compatible with that constraint

[TABLE]

where $g(x)$ is a positive, monotonic function that has to satisfy $g(x_{I}):=\sum_{i_{0}\in I}g(x_{i_{0}})$. For simplicity, we take $g$ as the identity, namely $g(y):=y$, which yields to .

[TABLE]

By applying the exponential to the above equation, it yields

[TABLE]

that yields the (off-diagonal) scale-invariant probability in Equation 12.

For the self-loops part, the steps are similar to the previous ones, except for the factor of $\frac{1}{2}$ that is introduced to avoid double counting. Indeed, referring to the network represented in Figure 9, Equation 36 yields

[TABLE]

which is solved only because the authors of [18] imposed the $\frac{1}{2}$ in the self-loop part.

Given the results we have just obtain, the defining relationship Equation 33 reads

[TABLE]

where the $p$ functional form is reported in Equation 12.

B.2 Parameter Estimation

In order to tackle the reconstruction problem -inferring the link probability starting from aggregated inputs-, the MSM parameters can be estimated by enforcing the degree sequence on average as done for the CM in Appendix A. For this reason, the output model is called the degree-corrected MSM (degcMSM). However, as the MSM is not an ERG model, the MoM and the maximum-likelihood estimation of the parameters may differ. For this reason, we will also consider the maximum-likelihood MSM (maxlMSM) to evaluate if, knowing the full structure of network, the latter procedure leads to a better modelling of the underlying structure. To stress the point, the comparison among the degcMSM and maxlMSM is a theoretical exercise the former assumes the availability of the degree sequence at best, whereas the maxlMSM needs the whole adjacency matrix, which implies the knowledge also of the degree sequence, but it is often unavailable due to policy restrictions [44].

B.2.1 Maximum Likelihood MSM

Assuming to know the adjacency matrix, the likelihood for the MSM reads

[TABLE]

Its gradient can be calculated by deriving the Equation 37 with respect to $x_{i}$ and $w_{i}$, namely

[TABLE]

which vanishes when the parameters are at the optimal point $\vec{x}^{*},\vec{w}^{*}$, i.e.

[TABLE]

where $p^{*}_{ij}:=p(\vec{x}^{*}_{i}\vec{x}^{*}_{j})\textnormal{ and }p^{*}_{ij}:=p(\vec{x}^{*}_{i},\vec{w}^{*}_{i})$. In particular, the second equation suggests that $\left\{w_{i}^{*}\right\}_{i\in[1,N]}$ are a function of $\vec{x}$ at the stationarity point, namely

[TABLE]

In addition, one can reduce the number of parameters by inserting the second equation into the first one, which leads to

[TABLE]

and, by integrating over $x_{i}$, to Equation 14.

In conclusion, it is enough to find the maximum of the off-diagonal likelihood Equation 14, and, then, fix the $\left\{w_{i}\right\}_{i\in[1,N]}$. This implies that the diagonal part of Equation 37 ($i=j$) vanishes, whereas the off-diagonal part is already maximized by means of the optimization procedure. Note that for a fair comparison among the CM (no self-loops modelling) and the MSM, the evaluation metrics (DEG, ANND, CC, $\dots$) don’t take into account the self-interactions.

The hessian was not calculated since it was used neither for a theoretical analysis nor in the numerical optimization -the SciPy methods can approximate it in a sophisticated way.

B.2.2 No natural statistical equivalence from structural equivalence

In the subsubsection III.3.1, we have seen that the degcMSM assigns the same $x$’s to nodes with the same degree. Here, since the likelihood sums over the links, one would expect that nodes with the same neighbors (structural equivalence - StructE) should have the same parameters (statistical Equivalence - StatE) at the optimal point Equation 42. However, this would not be the case (see next section), and one would need to enforce the equality among the parameters $x_{i}$ ultimately leading to the statistical equivalence of the StructE nodes. Roughly, by imposing the StatE, one discards a possible source of redundancy and lets the model focus on the backbone of the observed network, i.e. the inequivalent classes.

To show that the stationarity condition is not sufficient to guarantee the StatE, let us assume that the node $i,s$ shares the same neighbors, namely $N(i)\equiv N(s)$ where $N(i):=\left\{j:j\in[1,N],a_{ij}=1,j\neq i\right\}$ (self-loop included). This implies that their adjacency entries are the same, i.e. $a_{ij}\equiv a_{sj}\,\forall j$ as nodes $1,3$ in the Figure 10. Furthermore, from the stationarity condition, one can obtain the following relationships among $x^{*}_{i},x^{*}_{s}$, formally

[TABLE]

where we highlighted the parameters’ dependence of the model. The previous equations show that, even thought the nodes are structurally equivalent, their statistical equivalence is not guaranteed from the previous equations. Therefore, the StatE has to be imposed by hand for all the StructE nodes, namely

[TABLE]

where $S(e):=\left\{j:N(j)\equiv N(e)\right\}$.

As final remarks, note that the degcMSM assign the same parameter $x_{i}$ (no self-loop parameter is involved) for nodes with the same degree. This equivalence class is larger than the StructE used so far, and so the degcMSM will have a smaller set of available parameters than the maximum-likelihood ones. However, this is not the case for the CM since the degree sequence is its sufficient statistics, i.e. the maximum of the likelihood is recovered by imposing the degree sequence on average. Hence, the cardinality of the set of parameters is the same of fitting the CM via maximum likelihood or fitting the node degrees.

B.3 Removal of deterministic nodes

Network reconstruction assumes that the observed $\mathbf{A}$ is a realization of a Bernoullian random process. However, it may happen that some nodes are deterministic, i.e. fully-connected (FC) or disconnected (D) to the other nodes. Therefore, their behavior is trivially recovered by setting the hidden variable of the FC nodes to infinity or the D nodes to zeros (see Equation 12). Their optimized parameters would be $x_{FC}\to\infty$ and $x_{D}=0$, respectively implying that $\left\{p_{FC,j}\equiv 1\right\}_{j\in[1,N]}$ and $\left\{p_{FC,j}\equiv 0\right\}_{j\in[1,N]}$. In turn, $\textnormal{ Var}(a_{FC,j})=\textnormal{ Var}(a_{D,j})=0\quad\,\forall j$ implying that their sampled edges are fixed, and not contributing to the ensemble fluctuation. Concretely, we have hard coded their variables, as seen before, to account for their roles in the graph.

A similar reasoning applies to the nodes connected only to FC, the so called only2fc nodes ($\nu$). In particular, in the degree-corrected case (both Configuration Model and degree-corrected MSM), by assuming there is only one FC node and focusing on the FC/only2fc equations,

[TABLE]

which solution reads $p_{FC,j}\equiv 1\,\forall j(\neq FC),p_{\nu j}\equiv 0\,\forall j(\neq\nu,FC)$. Similarly, in the maximum likelihood setting, the likelihood of only2fc reads

[TABLE]

where $\ln\left(p_{\nu,FC}\right)=0$ for the unconstrained optimization of $\mathcal{L}_{FC}$.

In both cases, the only2fc parameters are undetermined. Setting for convenience $x_{FC}=\tau^{2}$ with $\tau\to\infty$, they require that

[TABLE]

which are trivially solved by setting $x_{o}\stackrel{{\scriptstyle!}}{{=}}1/\tau$ as the $x_{j}\approx\textnormal{ const}$ with respect to $x_{o}$. As the FC-parameters are degenerate, we decided to fix them after the optimization. In addition, this results in a more stable maximization of the likelihood, since one is manually setting the values for the trivial (disconnected entities) or degenerate nodes (FC and only2fc).

Lastly, this procedure doesn’t spoil the * renormalization rule*. Including an FC node to a block-node produces an FC block-node; whereas a vanishing parameter (disconnected or only2fc) gives the freedom to the other terms in the summation (see Equation 20) to determine the group parameter.

B.4 Algorithmic Complexity

To describe the higher levels, one needs a lower computational complexity by means of a scale-invariant model than a single-scale (SSM) one. Specifically, the computational complexity of the RHS of Equation 36 (the only one possible to coarse-grain an SSM) reads

[TABLE]

where $c$ is the complexity for the calculation of the connection probability $p$ for the referred model, e.g. CM, and $N_{0}$ are the number of structural inequivalent nodes444Due to this, $N_{0}$ will be model dependent. at the microscopic level $\ell=0$. In contrast, assuming that the “sum” operation is $\mathcal{O}(1)$, the complexity of the LHS reads

[TABLE]

since there are $N_{\ell}$ summation to get $x_{I}$555To account for the summation over the $w_{I}$, one should add a factor of $N_{\ell}$ and all the $\binom{N_{\ell}}{2}$ pairs at the level $\ell$. Hence, the ratio reads

[TABLE]

For example at level $\ell=3$ for the ING network and maxlMSM, $N_{3}=67,N_{0}=972$,

[TABLE]

giving an improvement of 3 order of magnitude using the summed MSM rather that the full evaluation of the coarse-graining probability.

B.5 Numerical Comparison between the CM and MSM functional form

Here, we evaluated the CM and MSM functional forms for $x\in[0,10]$. The crucial point is that, even if the CM is not theoretically scale-invariant (see subsection A.1), it may be numerically similar to the MSM. Actually, for small $x$, $p_{CM}\approx p_{MSM}\approx x$. Hence, the CM expected quantities may be similar to the MSM ones, spoiling the idea the CM is a single-scale model, which should be generalizable.

Appendix C Scores

C.1 Comparison Between Probabilities

For every level $\ell$, both CM and MSM give rise to 3 probability functions: the fitted $\mathbf{P}^{(\ell)}_{fit}$, summed $\mathbf{P}_{sum}^{(\ell)}$ and the coarse-grained $\mathbf{P}_{cog}^{(\ell)}$. Technically, the first twos have the same functional form, while depending on the fitted and the summed parameters, respectively. Yet, $\mathbf{P}_{cog}^{(\ell)}$ has a different connection probability for the CM-like models (subsection A.1), whereas it is equal to the summed probability for the MSMs (subsection B.1). To display these differences, we cross compared them to understand their hallmarks.

Firstly, in Figure 2, one can find how much $\mathbf{P}_{sum}^{(\ell)}$ differs from $\mathbf{P}_{cog}^{(\ell)}$ both for the CM-like and MSM models. Secondly, the comparison among the summed $\mathbf{P}_{sum}^{(\ell)}$ against fitted $\mathbf{P}_{fit}^{(\ell)}$ probabilities is shown in the insets of 7a. In particular, we plotted the $2D$ histogram of the density $\rho_{bin}$ of points inside each bin (the total number of bins is $30^{2}$), i.e.

[TABLE]

where $n_{\textnormal{bin}}$ is the total number of points inside that bin and $N_{\ell}(N_{\ell}-1)$ the number of pairs. Note that we colored the bins according to $\rho_{bin}$ (creating a heatmap) - the bigger the value, the brighter the color. Lastly, we didn’t evaluate $\mathbf{P}_{fit}^{(\ell)}$ against $\mathbf{P}_{cog}^{(\ell)}$ since the coarse-grained probability is functionally different from the others for the CM-like models (see subsubsection III.4.2) whereas $\mathbf{P}_{cog}^{(\ell)}=\mathbf{P}_{sum}^{(\ell)}$ for the MSM. Therefore, we used $\mathbf{P}_{cog}^{(\ell)}$ only in the first comparison.

C.2 Degree, Average-Nearest Neighbor Degree and Clustering Coefficient

The fundamental topological properties of a network are the degree, the average nearest neighbor degree (ANND) and the clustering coefficient (CC) [28]. Formally, each of such measurements is a function $Y(\mathbf{A})$ of an $N\times N$ adjacency matrix representing a graph $\mathbf{G}$.

Here, we compute them both in the observed network $\mathbf{A}$ and as expected by the model $\mathbf{P}(\mathbf{A}|\textbf{X})$. More precisely, the degree counts the number of edges that are incident to a node $i$, i.e.

[TABLE]

and its expected value is given by

[TABLE]

where $\langle\cdot\rangle$ denotes the expected value over the ensemble of graphs sampled from $\mathbf{P}(\mathbf{A}|\textbf{X})$. Moving two-hops away from $i$, the ANND reports the average degree of the neighbors of the node $i$, i.e.

[TABLE]

whereas its expected value reads

[TABLE]

where in the second passage we took advantage on the first order approximation $\mathbb{E}[\frac{X}{Y}]\approx\frac{\mathbb{E}[X]}{\mathbb{E}[Y]}$ (delta approximation) [36]. Lastly, the CC is defined as the ratio among the number of triangles of node $i$ and its number of wedges, namely

[TABLE]

whereas the expected one is

[TABLE]

C.3 Triangle Density

Inspired by [11], we computed the expected number of triangles for every model at disposal666In this essay, our objective was to model probabilistically the observed network rather than describing it exactly, namely the limit where $p_{ij}\equiv a_{ij}\,\forall i>j$.. Specifically, the expected density of triangles at a certain level $\ell$ is defined as

[TABLE]

where $\triangle(\mathbf{G}^{(\ell)}_{k_{i_{\ell}}\geq c})$ is the number of observed triangles (see Equation 56) calculated on the subgraph $\mathbf{G}^{(\ell)}_{k_{i_{\ell}}\geq c}$ composed by the nodes $I_{c}:={\left\{i_{\ell}:k_{i_{\ell}}\geq c\right\}}$ with degree lower (or equal) than a threshold $c$ [10, 11]. Its expected value reads

[TABLE]

where $i\in I_{c},j\in I_{c},k\in I_{c}$ and the probabilities $\tilde{p}_{ij}$ refers to the summed model $\mathbf{P}_{sum}^{(\ell)}$.

C.4 On the ensemble estimation

For more complicated measurements, e.g. the variance of the ANND, the delta approximation won’t be valid and one has to estimate them as the average over a sufficiently large ensemble $\mathcal{A}:=\left\{\mathbf{A}_{s}\right\}_{s\in[1,\mathcal{S}]}$ where $\mathcal{S}$ is the number of graphs. In particular, having optimized the parameters of the model, we can generate unbiased realizations $\mathcal{A}$ by sampling each $a_{ij}$ independently with probability $p_{ij}$ [36, 18].

In the limit of $\mathcal{S}\to\infty$, the sampled average of any measure $Y_{i}$ meets its analytical estimations $\langle Y_{i}\rangle$ [18], i.e.

[TABLE]

where $\mathbf{B}\in\mathcal{A}_{N}$ is a matrix drawn from the set of the undirected binary graphs $\mathcal{A}_{N}$ of $N$ nodes.

Lastly, to estimate the uncertainty of the model over the sampled realizations, we calculated the dispersion intervals $\Delta_{c}(\langle Y_{i}\rangle):=\langle Y_{i}\rangle\pm 2\sigma_{\langle Y_{i}\rangle}$ [42]. In other words, an interval centered in $\langle Y_{i}\rangle$, but $4\sigma_{\langle Y_{i}\rangle}$ wide.

C.5 Reconstruction Accuracy

In order to have a cross-comparison among all the levels and models, we exploited the reconstruction accuracy [42]. This measure is defined as the fraction of times an observed statistics $Y_{i}$ falls within the dispersion interval $\Delta_{c}(\langle Y_{i}\rangle)$ (see subsection C.2). More formally, the reconstruction accuracy at level $\ell$ for the statistics $Y$ is defined as

[TABLE]

where $\mathbb{I}$ is the indicator function777Further refinements are possible, but we stick with this definition for the sake of simplicity.. Roughly, it counts the frequency at which the sampled ensemble includes the observed statistics. If all the observed statistics, e.g. degrees, were included in the interval, the accuracy would be $1$, whereas the accuracy would be [math] if none of them were included.

C.6 Rescaled ROC and PR Curves

By considering the introduced models as binary classifiers it seems natural to inspect their expected confusion matrix, Receiver Operating Characteristic (ROC) and the Precision-Recall (PR) curves [45, 46]. In particular, the expected confusion matrix is a $2\times 2$ matrix that reports the expected value of True Positives (TP), i.e. $\langle TP\rangle:=\sum_{i<j}a_{ij}p_{ij}$, False Positives (FP), i.e. $\langle FP\rangle:=\sum_{i<j}(1-a_{ij})p_{ij}$, True Negatives (TN), i.e. $\langle TN\rangle:=\sum_{i<j}\left(1-a_{ij}\right)\left((1-p_{ij})\right)$, and False Negatives (FN), i.e. $\langle FN\rangle:=\sum_{i<j}a_{ij}\left(1-p_{ij}\right)$ [47]. By combining these scores, one recovers the True Positive Rate ($TPR$), the False Positive Rate ($FPR$) and the Positive Predictive Value ($PPV$) [46], namely

[TABLE]

where $Pos:=\sum_{i<j}a_{ij}=L$, $Neg:=\sum_{i<j}(1-a_{ij})=\binom{N}{2}-L$.

Instead of taking these expected values, another common practice is to set the entries s.t. $p_{ij}\leq\epsilon$ for a fixed threshold $\epsilon$, and otherwise to [math]. For convenience, we will call them $TPR(\epsilon),FPR(\epsilon),PPV(\epsilon)$. Furthermore, the ROC and PR curves are obtained by spanning $\epsilon\in[0,1]$ [45] and plotting the TPR-VS-FPR, and the TPR-VS-PPV, respectively (see 7b). Note that if $\epsilon\to 1$, then $TPR(\epsilon\to\infty)=FPR(\epsilon\to\infty)=0,PPV(\epsilon\to\infty):=1$, whereas if $\epsilon=0$, $TPR(\epsilon=0)=FPR(\epsilon=0)=1,PPV(\epsilon=0)=p:=\frac{Pos}{Pos+Neg}=\frac{L}{\binom{N}{2}}$.

Since the last point of the $PPV$ is the density, the “naive classifier” (NC) - always predicting the majority class - can have a high AUC-PR (On the ROC curve the NC always lies on the diagonal). For this reason, one can rescale the Area Under the Curve (AUC) of each model with respect to the AUC predicted by the NC. Specifically, one defines

[TABLE]

Therefore, the perfect classifier still have $\textnormal{ AUC-ROC }=\textnormal{ AUC-PR }=1$ but the random one $\textnormal{ AUC-ROC }=\textnormal{ AUC-PR }=0$. The new AUCs can be negative, as a signal of a worse performance than the random classifier. Furthermore, this novelty does not spoil the ranking of the summed models which was the ultimate objective of the AUC scores.

Appendix D World Trade Web results

In this section, we report the same results presented in the main text, but for the World Trade Web. The figure analysis are reported in the caption of each figure.