Identifying maximal sets of significantly interacting nodes in higher-order networks

We introduce a method for the detection of Statistically Validated Simplices in higher-order networks. Statistically validated simplices represent the maximal sets of nodes of any size that consistently interact collectively and do not include co-interacting nodes that appears only occasionally. Using properly designed higher-order benchmarks, we show that our approach is highly effective in systems where the maximal sets are likely to be diluted into interactions of larger sizes that include occasional participants. By applying our method to two real world datasets, we also show how it allows to detect simplices whose nodes are characterized by significant levels of similarity, providing new insights on the generative processes of real world higher-order networks.

Over the last decades we have witnessed an explosion in the number and size of data describing the different domains of the world we live in. In particular networks, collection of nodes and edges linking pairs of them, have emerged as a very valuable tool to describe and model the topology of all such real-world systems where interactions among agents play a key role, from trade economics to complex biological processes [1]. Nevertheless, mounting evidence is now suggesting that networks can only provide a limited description of many complex systems, as they are inherently limited to model pairwise interactions [2,3]. Higher-order interactions have been observed in a variety of contexts, from human interactions in face-to-face systems [4] to collaboration networks [5], ecological systems [6] and the human brain [7,8]. Interestingly, such non-dyadic interactions were shown to give rise to new collective phenomena in dynamics such as contagion [9,10] and cooperative process [11], as well as diffusion [12,13] and synchronization [14][15][16].
To properly encode such higher-order systems [17], more complicated mathematical frameworks are needed, such as simplicial complexes and hypergraphs. In hypergraphs, generalised hyperedges encode interactions among an arbitrary number of elementary units. In simplicial complexes higher-order interactions are described by simplices that satisfy a particular inclusion property, such that also all subfaces of a given simplex must exist. Practically, this means that if n interacting nodes are described by a simplex, also all interactions among all possible subsets of nodes are considered to exist. This additional constrain comes at the benefit of enabling the use of computational techniques from algebraic topology, typically known as topological data analysis, to efficiently describe the topological features of such systems [18,19]. Given the growing interest in higher-order interactions, over the years several network metrics have been generalised to account for the presence of non-dyadic interactions. These include local scores such as centrality measures [20][21][22], higher-order motif analysis [23,24] and simplicial closure [25], the proposal of methods to extract higher-order community structure [26][27][28][29], or unveiling higher-order dependencies in multivariate time-series [30].
The rising interest towards higher-order complex systems calls for the design of suitable filtering techniques. In fact, in light of the explosion of data on systems belonging to different domains, filtering methods have attracted attention as tools that allow to focus on relevant connections that are highly representative of a system, while discarding all the redundant and/or noisy information. This approach is useful to detect relevant trends in the data, otherwise hidden by noise or by the intrinsic heterogeneity of investigated complex systems.
Filtering approaches for complex systems can be divided into two main families. The first family of filtering approaches relies on imposing topological constraints that force to reduce the number of connections in the filtered system. The first example of such approaches is the Minimum Spanning Tree, which reduce a system of n connected nodes to a tree with n − 1 links that are placed starting from those with the highest weight [31]. This approach has been extended by proposing new topological constraints [32,33] and the robustness of the filtering procedures has been assessed through statistical procedures [34,35].
In the second family we find methods that drop the connections that are compatible with an appropriate null model. Null models are used to produce ensembles of replicas of the system, where some key feature are preserved while the generative model that produce links and/or hyperlinks is randomized. An example is the configuration model, in which the degree of each node is exactly preserved but the edges are placed randomly. Different null models have been proposed in the last decade, for different classes of systems. There are null models for weighted networks [36][37][38], bipartite networks [39][40][41][42][43][44], temporal networks [45] and correlation matrices [46,47]. Recently, the first null models for higher-order systems have been proposed [45,48].
Statistically validated filtering of complex systems has been applied to several domains, including actor/movies networks [39], trading decisions of investors [49][50][51][52], financial transactions between major market agents [53], mobile phone calls of a large set of users [40,54], financial credit transactions occurring in an Interbank market [41], intraday lead-lag relationships of returns of financial assets traded in major financial markets [55], the Japanese credit market [56], features of communities detected in complex networks [57],international trade networks [58], social networks of news consumption [59], and rating networks of e-commerce platforms [60].
In this work, we propose a filtering methodology for higher-order complex systems that is able of extracting Statistically Validated Simplices from collectively interacting groups of nodes. Specifically, we develop a method using a null model designed for repeated higher-order interactions, that is able to identify the maximal sets of nodes that consistently interact collectively and discard the nodes that co-interact with them only occasionally. We label these minimal sets of nodes as Statistically Validated Simplices, as the filtered significant units share the inclusion property that defines a simplex. Since our null models expands a class of existing models [39,48], it inherits from them the capability of preserving the heterogeneity of the system. The paper is organized as follows: (i) in section I we introduce and motivate our method and discuss a possible approximation of the proposed p-value; (ii) in section II we introduce a class of synthetic benchmarks that we use to test the accuracy and precision of our method against available alternative in literature. We also apply our method to empirical data and show that it is able of extracting relevant information on the investigated systems; (iii) finally in section III we draw our conclusions.

I. FILTERING METHOD
We aim at obtaining a tool that is capable of identifying groups of agents that are consistently interacting among them. Specifically, we want our groups to contain agents that are consistently appearing together in collective interactions, without including those agents that might appear in some of these interactions but whose presence is occasional. For example, consider a set of n products that are consistently bought together in the same basket because of functional relatedness. Given the heterogeneity of shopping habits, it is likely that these baskets will feature occasional other products. We want our tool to identify the n products that constitute a cohesive unit and discard all those other ones that are copurchased only occasionally. In the following we adopt the formalism of higher-order networks [2], that naturally provides the concepts required to build our tool. Within this formalism, agents are mapped into nodes and higherorder interactions between them are mapped into hyperlinks. We also recall the concept of simplicial complex as a collection of simplices, where a simplex of order n is a set of n + 1 interacting nodes, with the condition that all interactions between any subset of the nodes also exist. In what follows, to avoid confusion, we will always characterize hyperlinks and/or simplices through their size, which is the actual number of nodes that are interacting. A hyperlink/simplex of size n corresponds to an interaction between n nodes. Thus in this framework our research question can be rephrased into the problem of transforming a collection of hyperlinks that encode recorded higher-order interactions between a set of real world agents into a simplicial complex where each simplex represent a set of statistically validated interacting nodes. This transformation comes with significant implications. Indeed, a simplicial complex is a mathematical tool that allows to analytically generalize to higherorder systems many results that hold for pairwise networks [9,61]. However, not all real world collections of hyperlinks constitute a simplicial complex, as a simplicial complex K is valid only if, for every simplex of size n that constitutes it, all its sub-simplices of any size n < n also belong to K. A collection of real world hyperlinks in general does not hold this property, as there is no guarantee that for each set of n nodes interacting in a hyperlink all subsets of n < n nodes will also be interacting. On the contrary, we want our statistically validated groups of co-interacting nodes to be simplices because we consider all subgroups of n nodes that are significantly interacting together to be also significantly interacting. To keep our example, if n products are consistently bought together, any subgroup of these products will also be significant, as they will appear on the same baskets. In other words, we want to identify the statistically significant higher-order behavior of agents -i.e. the core of n products that are bought together consistently -, going beyond the detection of higher-order observables -i.e. all the specific combinations of core products and occasional ones that are sometimes bought together without constituting a significant basket (Fig. 1). In what follows, we will call the filtered set of significant sets of interacting nodes Statistically Validated Simplices (SVS). Given n nodes that are active in a system, we need a statistical criterion to assess whether the number of times in which they interact together is statistically significant. The simplest way to obtain this result would be to set a threshold on the number of times required to consider a group significant. However, this approach is subopti-mal because it is incapable of taking into account the heterogeneity in activity among agents. In fact, as the individual activity of agents in complex systems is usually strongly heterogeneous, a single threshold -e.g. one extracted from the distribution of activity of all nodeswould be too small for hyperlinks that contains the most active nodes and too large for hyperlinks that contains the least active nodes. Also, it would require to design an arbitrary set of rules to deal with hyperlinks whose participating nodes have heterogeneous levels of activity. For this reason, the validation of networks and hypergraphs is usually performed by means of a null hypothesis that is able to naturally incorporate the heterogeneous activity of nodes. Here we introduce a modified version of the null hypothesis that was designed for validating hyperlinks for Statistically Validated Hypergraphs (SVH) [48]. In fact, the SVH approach is limited to validate higher-order observables of a specific size n, i.e. observed hyperlinks of n nodes. The statistical significance of a hyperlink is assessed by considering its weight X -i.e. the number of times in which it is observed -, the number of hyperlinks of size n to which each node participates and the global number of hyperlinks of size n. From these parameters one can extract a p-value that represents the probability of observing the same weight X or more under a null hypothesis in which each node randomly selects its n − 1 partners. Here, we want to modify this approach by considering and testing all possible combinations of nodes that co-interact -no matter with how many other nodes -at least once. To obtain this, we modify the null hypothesis used in the SVH in the following way.
For simplicity, we start from the case of 3 nodes i, j, k that co-interact N ijk times in hyperlinks of size n >= 3. The three nodes appear respectively in N i , N j , N k hyperlinks. As we are not limiting anymore the possible hyperlinks to those whose size is exactly 3, we now set as N the total number of hyperlinks observed in the system. This choice for the value of N consider all hyperlinks equally likely to contain any number of nodes. Alternatively, N can be set as the number of hyperlinks with size larger or equal n -in our case n = 3. The hyperlinks that are counted in N ijk all contains the nodes i, j, k but, differently from the case of SVH, are not limited to them. For example, a hyperlink that contains the nodes i, j, k, l also contributes to N ijk . Under the null hypothesis that each node selects randomly the hyperlinks to which it participates -and thus its n − 1 counterparts in any hyperlinkthe probability of observing i, j, k interacting N ijk times is where H(N AB |N, N A , N B ) is the hypergeometric distribution that computes the probability of having an intersection of size N AB between two sets A and B of size N A and N B given N total elements. The probability p(N ijk ) in Eq. 1 represents the probability of having a random intersection of size N ijk between the three sets of hyperlinks of nodes i, j, k out of N total hyperlinks [62], and is obtained through the convolution of two instances of the hypergeometric distribution. Indeed, we start from the probability of having an intersection of size X between nodes i and j and multiply it with the probability of having an intersection of size N ijk between the hyperlinks of node k and the intersection set of size X between the hyperlinks of node i and those of node j. This product is then summed over all possible values of X, i.e., all possible intersections between i and j which are compatible with the observed number of interactions between all the three nodes, which lie in the interval Eq. 1 is symmetrical around permutations of i, j, k as the specific order with which they are considered does not affect the value of the probability. Starting from Eq. 1, we then compute a p-value for the symplex that contains i,j and k through the survival function, The p-value provides the probability of observing under our null hypothesis N ijk or more hyperlinks that contain -but are not limited to -the nodes i, j, k. For a generic simplex of n nodes, Eq. 1 becomes ... × In order to assess the statistical significance of a simplex of size n, its p-value is tested against a threshold of statistical significance α, after adding a multiple hypothesis test correction which is needed because of the high number of tests -one per each simplex. In all the results presented in this paper we use α = 0.01. However, since we consider all smaller combinations of nodes constituting a significant simplex to be themselves significant simplices, we proceed to p-value testing in a sorted order. In fact, we start from testing the largest simplex and we then proceed towards the smallest. If a simplex of size n passes our statistical test and is thus selected as informative because it rejects the null hypothesis, we do not test any of its smaller subsets. Moreover, as the number of possible simplices of size n from N nodes nodes explodes as N nodes n , we only test simplices whose n nodes appear in at least one observed hyperlinks. This means that we consider candidates for being significant simplices only sets of nodes that interact at least once. Thus, for each size n, we proceed from largest to smallest (where the largest n corresponds to the size of the largest observed hyperlink), and we test all candidates size by size after correcting the significance threshold by means of the control for the False Discovery Rate [63]. The total number of test considered at each size n is N t = N nodes n which is the number of all possible combinations with no repetitions of size n of the N nodes elements that are active in the system. Thus, when applying the control for the False Discovery Rate method we start from a Bonferroni threshold computed as α B = α/N t .

A. Approximated p-values
Extracting the probability of Eq. 3 requires the convolution of n − 1 instances of hypergeometric distributions, each of which is built through the combination of several binomial coefficients. For this reason, this task becomes computationally intense for larger n. On top of this, our approach may require the computation of several p-values over the same nodes, as any set of n nodes which does not represent a significant simplex is split into n candidate simplices of size n − 1. If not significant, each of these complexes is further split into candidate simplices of size n − 2, and so on, quickly causing the number of computationally intense p-value that need to be computed to explode. To overcome this computational bottleneck, we have tested an approximated version of the p-value. Eq. 3 represents the probability that n nodes appear in N 12...n hyperlinks after each of them has randomly sampled its hyperlinks from an urn. As the random sample performed by each node is independent from the others, we can approximate this estimation using the binomial distribution, turning Eq. 3 into Ni N . In other words, we are estimating the probability of observing N 12...n hyperlinks (out of the N total hyperlinks observed) that involve the n nodes, where each hyperlink has a probability of being observed which is the product of the normalized number of hyperlinks Ni N in which each node appears. Our approximated p-value recalls the one proposed in Ref. [45], as both approaches use the binomial distribution. However, here we model differently the probability of observing a collective interaction. In Ref. [45], given the specific nature of the investigated data that only tracks temporal pairwise interactions, the p in the binomial distribution is obtained by multiplying the probabilities of observing all pairwise links that constitute a clique of the same size of the tested hyperlink. Here, we obtain it by multiplying the normalized number of hyperlinks of each node involved in the interaction. The approximation of Eq. 4 lies in the fact that we do not take into account anymore the n-dimensional array of individual activities of each node in the hyperlink {N i }, but we consider the normalized 1-D product between them.
We thus explored the precision of this approximation in different scenarios (Fig. 2). To do so, we randomly generated groups of nodes of different order with different values of {N i } (out of N total hyperlinks) and different values for the number of collective hyperlinks N 1,2,...,n . Specifically, for each size n we extracted n values of N i from the uniform distribution U (1, f N ), where f is the fraction of maximum number of hyperlinks out of the total N that we allow a node to appear in. We then extracted N 1,2,...,n from the uniform distribution U (1, min(N i )). We then used the values (N 1,2,...,n , N i , N ) to extract a p-value with the exact approach of Eq. 3 and the approximated one of Eq. 4. We varied the size of the simplices to validate in the interval [1,4]. Each column of Fig. 2 presents results with a different size. We also varied the total number of interactions  Fig. 2 presents results with a different value of f . We find that across the different settings, the approximated p-value tends to overestimate the effective probability that the number of interactions within a group is compatible with our null hypothesis. This can lead to False Negative, i.e. to miss groups that are significant but for which we overestimate the actual p-value. However, this overestimation is pronounced only for values of f close to 1, which means cases in which a system is made of a small amount of nodes each of which is responsible, on average, of a significant portion of the complete number of hyperlinks. In cases in which a single node is at maximum present in a smaller fraction of the overall interactions, the approximation works well and the p-values computed using the two approaches fall on the diagonal. In all the applications that we are going to showcase in the following, we have used the approximated p-value after verifying that the maximum value of f for each node in all investigated systems was below 0.1.

A. Definition and investigation of benchmarks
To test the effectiveness of our approach, we create benchmark hypergraphs made of N nodes in which we implant T simplices for each size n. In this work, we explored N = {100, 200, 500} and T = {100, 200, 400}, with n ranging in the interval [2,4]. For each simplex, we create l hyperlinks, where l is sampled from a Binomial distribution with p bin = 0.5 and n bin = 20. Each hyperlink contains the n nodes of the simplex plus a randomly selected set of the other N − n nodes whose exact size is sampled from the uniform distribution U (0, n max ). In our simulations, we set n max = 6. In this way the n nodes of a simplex interact l times in different hyperlinks, which always contain the n core nodes plus a set of randomly selected ones. Moreover, we introduce a pa- rameter c that controls the level of simplicial closure at each size. Specifically, c is the fraction of groups of size n that constitute a clique of size n + 1. For n = 2, c corresponds to the concept of triadic closure and represents the fraction of dyads that are closing a triangle. For each realization of our benchmark, we use (i) the pairwise method of Statistically Validated Networks (SVN) [39], (ii) the method of Statistically Validated Hypergraphs (SVH), that only validates observed hyperlinks of fixed size [48], and (iii) our new method of Statistically Validated Simplices (SVS) to extract the simplices rejecting the null hypothesis. To validate our results, we compute the True Positive Rate (TPR), that quantifies the ratio of true groups that we are able to identify, and the False Discovery Rate (FDR) that quantifies the ratio of erroneously detected groups that the different methods select. Specifically, TPR= T P T P +F N and F DR = F P F P +T P , where TP are the True Positives (in our case, simplices correctly identified), FN are the False Negatives (significant simplices present in the system but not identified by our method) and FP are the False Positives (simplices wrongly identified as significant). Fig. 3a shows TPR as a function of c for N=200 and T=100, but the observed behavior is independent of the specific values of N and T. Each point is the median of 100 benchmark realizations. The color bands represent the interval between the 10 t h and the 90 t h percentile. The performance of SVS and SVH is not affected by c, as both are methods designed to take into account the higher-order nature of interactions. In fact, if 3 nodes i, j, k constitute a clique of (significant) pairwise interactions but never interact altogether, SVS and SVH will not statistically validate the corresponding simplex because N ijk = 0 even if N ij , N ik and N jk are not. As a consequence, SVS and SVH will not miss {ij}, {jk} and {ik} as significant simplices. Conversely, the performance of SVN significantly decreases with c specifically because the SVN puts together nodes that never interacted altogether and misses the smaller true simplices. In our example, SVN will erroneously identify {ijk} as a significant simplex and will miss that the significant ones are the couples {ij}, {jk} and {ik}. We also observe that SVS significantly outperforms SVH. This happens because in our benchmark hypergraphs we allow the n nodes of a simplex to interact in larger groups. While SVS counts all the hyperlinks in which the n nodes appear, the SVH approach misses all those in which there are additional, randomly selected nodes. For this reason, SVH is more likely to underestimate the significance of a simplex and thus have a lower value of TPR. Fig. 3b plots the FDR as a function of c. SVH and SVS almost never detect wrong groups as they distinguish between interactions of different order. SVN instead has a significantly higher value of FDR, which worsens for high c because the number of simplices erroneously aggregated at larger sizes increases, as pointed out in our example discussed previously. To better characterize the difference in performance between SVS and SVH, we add to the benchmark a parameter f that regulates the fraction of times in which the n nodes that constitute a simplex interact in hyperlinks of size n > n. The larger f , the more we expect SVH to underestimate the number N 12...n of collective interactions and thus will likely fail to validate the true simplex. In Fig. 3c we show TPR as a function of f . Indeed, we see that SVH significantly worsens its performance when f increases. Fig. 3d shows the behavior of FDR as a function of f . SVH does not detect a significantly larger ratio of false simplices until f is so large that it detects very few positives (compare with TPR for f ≥ 0.8, all of which are false.
B. Application to real data: Baskets of co-purchased goods at Walmart We analyse a dataset that tracks 69,906 co-purchases of 88,860 products at Walmart [64]. This system can be mapped into a hypergraph where hyperlinks represent copurchases. In this context, we use the SVS approach to extract simplices of products that are inherently bought together. When extracting the Statistically Validated Simplices we exclude from the hypergraph all baskets of size larger of 10. Indeed, each large basket of size k which is not statistically significant produces a large number of candidate simplices of smaller size k -specifically it produces k k candidates. Given that larger baskets have small number of co-occurrence, we cut the system at size 10 in order to avoid the explosion in the number of hyperlinks with low chance of getting validated that we check. Fig. 4a shows the number of validated (blue) and not validated (orange) simplices for each size. The number of validated simplices exponentially decays with size, and, for large n, it is several orders of magnitude smaller than the simplices that are observed at least once but are not rejecting the null hypothesis. The dataset also reports the category of each product. By using this information, we define a proxy of "purity" of co-purchased goods. Purity is a measure that goes from 0 (if all products in a group belong to different categories) to 1 (if they all belong to the same category), and is defined as P = n n−1 (1 − C n ), where n is the size of the simplex and C is the number of unique categories to which the n products belong. Fig. 4b shows the "purity" with respect to category classifications of the validated (blue) and not validated (orange) simplices as a function of their size. In the figure, the solid lines plot the mean and the transparent bands the error of the mean. We find that the SVS of any size are more homogeneous with respect to category than the ones that are not validated.

C. Co-traded stocks belonging to the S&P100 index
We investigate a dataset tracking tick-by-tick transactions performed by 96 stocks that belong to the S&P100 index and are traded on any US venue in the month 12/2020 [65]. The data comes with a temporal resolution of 1 ms. By applying our method to this dataset, we detect SVSs of stocks that are co-traded in the investigated period in excess to what expected by a null model taking into account the heterogeneity observed in stock trading. Also here, before extracting the SVSs we exclude from the hypergraph all simplices larger than 10 stocks. We perform our analysis at a daily scale and then aggregate the results obtained for different days. Fig. 5a shows the number of validated (blue) and not validated (orange) simplices for each size. Both for validated and not validated simplices, Fig. 5b shows the mean purity with respect to the economic sector of stocks (solid lines) and its error (transparent bands) as a function of size. We find that validated simplices (blue line) have higher purity at smaller size but the difference disappear in larger ones. Fig. 5c shows the mean correlation (solid line) and its error (transparent bands) between couples of stocks in validated (blue) and not validated (orange) simplices. Correlation is computed on time series of financial returns extracted at a resolution of 1 minute. We find that stocks belonging to validated simplices have higher correlation patterns than those that are not validated.

III. DISCUSSION
The recent surge in interest towards higher-order networks, together with the increased availability of large amount of data on several real world systems, have lead to the introduction of the first filtering tools for hypergraphs [45,48]. However, the task of detecting statistically validated units of collectively interacting agents is multi-faceted, and the multidimensional nature of higherorder interactions calls for the design of different validation approaches. The existing methods have focused on the extension of validating approaches from pairwise in- teractions to hyperlinks [45] and the detection of statistically validated groupwise interactions that occur only at a specific size [48]. Here, we introduced a different and more general framework for the detection of Statistically Validated Simplices, which represent sets of nodes that appear together in interactions of any size. The computational efforts needed for the statistical validation are demanding but accessible with ordinary PCs for groups up to a size of about 10 elements. To speed up the computation we have proposed an approximated procedure that work well in a large number of cases as shown in our investigation of benchmarks. Using properly designed hypergraph benchmarks, we show that our approach is highly effective on higher-order systems (i) in which informative interactions among n core nodes are likely to be diluted into larger hyperlinks that feature occasional nodes and (ii) with a high degree of simplicial closure. We also show how, on real data, our approach allows to detect simplices whose nodes are characterized by high levels of similarity, which might help in shredding light on the generative process of the underlying hypergraphs.

Code Availability
The code to extract Statistically Validated Simplices is available at https://github.com/musci8/SVH.


extfieldsofstudy:
Computer Science
Engineering
Mathematics

text:
Distributed Constraint-Coupled Optimization over Lossy Networks

This paper considers distributed resource allocation and sum-preserving constrained optimization over lossy networks, where the links are unreliable and subject to packet drops. We define the conditions to ensure convergence under packet drops and link removal by focusing on two main properties of our allocation algorithm: (i) The weight-stochastic condition in typical consensus schemes is reduced to balanced weights, with no need for readjusting the weights to satisfy stochasticity. (ii) The algorithm does not require all-time connectivity but instead uniform connectivity over some non-overlapping finite time intervals. First, we prove that our algorithm provides primal-feasible allocation at every iteration step and converges under the conditions (i)-(ii) and some other mild conditions on the nonlinear iterative dynamics. These nonlinearities address possible practical constraints in real applications due to, for example, saturation or quantization among others. Then, using (i)-(ii) and the notion of bond-percolation theory, we relate the packet drop rate and the network percolation threshold to the (finite) number of iterations ensuring uniform connectivity and, thus, convergence towards the optimum value.

I. INTRODUCTION
Distributed algorithms are widely considered for optimization and learning over networks with several seminal works addressing the case of networks with reliable links [1]- [10]. Due to the unreliable nature of wireless communication, these results are not readily applicable in wireless networks. In such networks, randomness of the links should be considered as link-failures due to packet dropouts are inevitable. The unreliability of information exchange in such settings can influence the algorithm's convergence and solution feasibility.
Equality-constraint optimization problems intend to optimize the allocation cost subject to constant sum of overall resources. The existing literature either focus on dual-based methods and ADMM solutions [1]- [5] or primal-feasible gradient-tracking methods [6]- [10]. These allocation algorithms are developed over reliable communication networks and require either (i) stochastic weight design over the network or (ii) all-time network connectivity as a giant component [7], [8] or both [1]- [6]. However, in unreliable networks with packet loss and message drops over the communication links, these restrictive conditions do not necessarily hold. As the messages are not delivered over some links at different times, the network resembles a dynamic graph topology over time that may even lose its connectivity at some iterations, failing condition (ii) for high drop rates. On the other hand, even assuming a connectivity guarantee, condition (i) mandates robust consensus algorithms (as in [11], [12]) to modify the stochastic weights on the shared information after the drop of some messages.
The proposed allocation algorithm in this paper relaxes the stochastic condition in [1]- [6] on the adjacency matrix to being only a balanced weight matrix. Assuming the common knowledge of the delivered messages over the undirected links (i.e., both sender and receiver know whether its packet is delivered) [11], [13], there is no need for compensation strategies to readjust the weights to make them stochastic again in case of link removal or packet drops. Further, our solution only needs uniform-connectivity over time relaxing the all-time connectivity requirement over undirected [1]- [5], [7] and directed networks [6], [8]. This implies that the network can lose path-connectivity (i.e., connectivity over a sequence of directly-connected nodes) between some nodes at some iterations.
In this paper, the loss of packets at every iteration is modeled as removing the associated links over the network. Therefore, for high packet drop rates, a connected network may go through a transition phase and lose its connectivity (due to link failure). This is studied via percolation theory in network science literature [14] and is discussed in Section II-A in details. In particular, bond-percolation refers to phase-transition (or percolation) in the network connectivity under certain rate of link failure or removal [15]. Similar threshold-based approaches are adopted in the epidemic processes with a spread of a certain disease or virus over the network [16], [17]. As a typical approach, we analyze this bond-percolation analysis over general random networks, where the randomness stems from the link failures, e.g., due to random packet dropouts in wireless sensor networks [18]. Some well-known existing random models include Erdos-Renyi (ER), Scale-Free (SF), and Small-World (SW) networks. Such models are further known to be a typical representative of networks encountered in many real-world applications. Therefore, the literature focuses on studying such networks' properties to understand and resemble the behaviour of a diverse range of practical large-scale networks, including IoT, social, financial, and transportation networks. We relate the packet drop rate to the uniform connectivity of the network and, in turn, its bond-percolation threshold such that the convergence of the optimization algorithm is guaranteed.
Paper Organization: Section II formulates the problem with some background on the bond-percolation theory. Section III provides the proposed allocation algorithm and its convergence analysis. Section IV states the convergence under packet drops. Section V and VI provide simulations and concluding remarks.

A. Background on Graph Theory
In many networked applications, the topology and interactions of entities resemble a graph model. In such a network topology, connectivity plays a vital role in the convergence of the adopted algorithms. Any change in the network may cause a transition from connected to disconnected, which is studied via percolation theory.
Definition 1: [14] Given the probability of link removal equal to p, bond-percolation is defined as the probability threshold p c such that for p > p c there is no giant connected component in the network with probability 1, i.e., the network loses its connectivity. Recall that, following Kolmogorov's zero-one law, one can also claim that for probability of link removal p < p c the network preserves its connectivity [14]. In other words, this critical probability p c represents a two-sided phase transition point in terms of network connectivity.
In the perspective of packet drop analysis, this percolation threshold is similarly defined as the probability where there exists no reliable delivery-path between (at least) two nodes over the network, where a path between nodes i, j denotes a sequence of linked nodes starting at i and ending at j. Due to its complexity, no analytical solution exists to define the bond-percolation threshold for general networks, and the existing literature is mostly limited to lattice and grid networks. For example, some rigorous bounds for different lattice graphs are given in [14], [19]. However, many experimental and numerical works on this problem exist in the literature based on the Monte-Carlo simulation [15]. For a survey of percolation theory over complex networks and wireless networks see, for example, [19] and [20].
On the other hand, random graph models are analyzed to estimate the percolation thresholds of similar large-scale and complex real networks. An interesting result is reported following the nearest neighbour rule and the well-known scale-free (SF) model, assuming a set of nodes scattered via the Poisson point process, for example, in networks of cellular or mobile phone base stations. One can show that connectivity of nodes to (at least) m = 3 nearest neighbours, in the presence of lossy links and drop-outs, guarantees that the entire network remains connected [20]. For the ER random graph models [21], the bond-percolation threshold is defined equal to 1 N with N denoting the average node degree. For the Small-World (SW) networks, this threshold is exactly determined in [22]. Table I summarizes the bondpercolation thresholds for some typical networks.   TABLE I: Bond-percolation threshold for different networks   Network type Bond-percolation threshold p c Square Grid [15] 0.5 ER [21] 1 N with N as mean degree with θ as shortcut probability SF [19] with ζ(·, ·) as Hurwitz zeta function, N m as min node degree, and σ as power-law degree Next, we provide some relevant notions on algebraic graph theory [23], [24] and Laplacian analysis over graphs [25]. Denote by W the adjacency matrix associated with the network G = {V, E}, with weights W ij > 0 as the weight on the link (j, i) ∈ E and 0 for (j, i) / ∈ E. If the network is undirected with balanced weights, W is symmetric. For such a network, define For a connected network (containing a spanning-tree) with symmetric L, its eigen-spectrum includes one and only one 0 eigenvalue associated with the eigen-vector 1 n . Recall that spectral localization of the L matrix plays a key role on the convergence properties of the consensus algorithms and distributed optimization methods, see more details in [25].

B. The Equality-Constraint Allocation Problem
The constraint-coupled optimization problem in this work is defined as follows: Intuitively, solving P 1 gives the optimal allocation of resources x for which the cost F (x) is minimized. Recall that the equality-constraint is known as the feasibility constraint or the sum-preserving constraint.
In general, there might be some local box constraints m i ≤ x i ≤ M i involved in P 1 that can be addressed using additive penalty terms [26], [27] or the so-called barrier functions [2]. Proper initialization for x i (0) under such box constraints is discussed in [8]. Following the KKT conditions and under Assumption 1, it is clear that for the unique optimal point x * we have ∇F (x * ) = ϕ * 1 n with 1 n as the column vector of 1s and ∇F (x * ) := (∂f 1 (x * 1 ); · · · ; ∂f n (x * n )), assuming that this x * is in the range of the box constraints.

III. THE ALLOCATION ALGORITHM
We consider the following gradient-tracking dynamics to solve P 1 : with g n (·) and g l (·) representing some nonlinear functions on the nodes and the links. The following assumptions hold throughout the paper. Assumption 2: All the links over the network G(k) are undirected and weight-balanced, i.e., the adjacency matrix W (k) is symmetric. Further, there exists a B ∈ Z ≥0 such that the network G B (k) = k+B k G(k) is connected, i.e., the union network G B (k) includes a spanning-tree for all k ≥ 0 (implying uniform connectivity or B-connectivity).
Given the network structure, its associated balanced weight matrix W can be designed in a distributed manner using the strategy in [28]. In case of quantized updates, similar to allocation in [6], [29], integer weight-balancing strategy [30] can be used, however, with convergence to the ε-neighborhood of the optimizer [29].
Assumption 3: The nonlinear mappings g n , g l : R → R are odd and "strongly" sign-preserving with g(z)z > 0 for z = 0, g(0) = 0, and lim z→0 g(z) with κ l , K l and κ n , K n defined for g l (·) and g n (·), respectively. Such g n (·) and g l (·) include, for example, all monotonically increasing and Lipschitz functions. It can be proved that under Assumptions 2-3, the gradient-tracking solution (2) satisfies the feasibility constraint at all times k. This is called all-time feasibility, and, in resource allocation perspective, implies that the assigned resources n i=1 x i always meet the demand b and is a privilege of our gradient-tracking solution (2). Note that violating this constraint can damage or disrupt the service in many applications [2]. Our proposed allocation scheme is summarized in Algorithm ??.
Data: Input: W , N i , η, f i (·) Result: Output: Final state x(k) and cost F (x(k)) Initialization: Every node i sets k = 0 and randomly chooses m i ≤ x i (0) ≤ M i satisfying the feasibility while termination criteria NOT true; do Each node i receives g l (∂f j (x j (k)) from j ∈ N i ; Updates x i (k + 1) via Eq. (2) ; Shares ∂f i (x i (k + 1)) with its neighboring nodes j ∈ N i ; Sets k ← k + 1 ; end Algorithm 1: The Resource Allocation Algorithm Example applications in economic dispatch problem [9] and CPU scheduling and battery reservation over the smart grid [10] are given for g n (z) = z|z| v1−1 + z|z| v2−1 with 0 < v 1 < 1, v 2 > 1, and logarithmic quantizer q l (z) = sgn(z) exp(q u (log(|z|, ρ))) with q u (z) = ρ z ρ as the uniform quantizer, sgn(·) as the sign function, and [·] as rounding to the nearest integer. In case the function g n (·) or g l (·) is sign-preserving, but violates condition (3), one can guarantee convergence to the ε-neighborhood of the optimizer x * ; for example, with uniform quantization [31] and single-bit data exchange scenarios [32].
Theorem 1: Let Assumptions 1-3 hold. The dynamics converges to the optimal solution of P 1 for η(B + 1) < η := κnκ l λ2 uλ 2 n K 2 n K 2 l with λ 2 , λ n as the smallest non-zero and largest eigenvalue of G B (k) for all k ≥ 0.
Proof: We provide the sketch of the proof here. Uniqueness: Following Assumption 1 and the KKT conditions one can prove that there exists a unique optimizer x * to P 1 satisfying ∇F ∈ span(1 n ) [7].
Following the strong-convexity in Assumption 1 we have [26] F (x(k + B)) ≤ F (x(k)) We first consider the case B = 1. To prove F (k+1) ≤ F (k), we need to show where we dropped the k for notation simplicity. Recall the definition of the Laplacian-gradient tracking dynamics in [8] and define dispersion parameter ξ(k) := ∇F (k)−1 ⊤ n ∇F (k). Using the results in [25, Section V], from (2) and some mathematical manipulations based on Assumptions 2-3 it is sufficient that with the strict inequality for and for ξ = 0 n holds the equality implying ∇F ∈ span(1 n ). For B ≥ 1, the right-hand-side of (8) changes to η(B + 1) with parameters λ 2 , λ n defined for G B . This completes the proof. Note that, in the case of versatile network topology, one can choose the ratio λ2 λ 2 n for the spanning tree contained by G B (k) as the minimum connectivity ensured by Assumption 2. This min spanning tree problem could be a promising direction of future research to determine the bound on the convergence rate of the algorithm.
Some relevant results on the eigen-spectrum of L are given in the literature, for example, to estimate the spectral range based on the node degrees [25], its variation under link removal/addition [17], bounds on the algebraic connectivity λ 2 (G) ≥ 1 ndg with d g as the network diameter [24, p. 571]. The following remark distinguishes this work in terms of packet-drop tolerance.
Remark 1: In the proposed allocation dynamics (2): 1) In contrast to [1]- [7], we do not require the weight matrix to be stochastic, but only to be symmetric. Therefore, there is no need for weight compensation strategies [11], [12] after packet drops or any change in the network. 2) The algorithm converges under uniform connectivity, in contrast to all-time connectivity requirement, e.g., in [6]- [8] or ADMM-based solutions [1]- [5]. Therefore, although the network may lose its connectivity over some time periods, due to a high rate of packet drops or switching over sparse topologies, uniform connectivity over every B ∈ Z ≥0 steps is sufficient for convergence. This remark motivates the application over lossy networks, as discussed next.

IV. PACKET DROPS AND SPARSE CONNECTIVITY
In this section, we focus on networks with unreliable links. This can represent various scenarios, for instance, when nodes are activated from sleep mode at a random times as is common in devices with energy harvesting capabilities [33]. In general, data transmission over wireless networks is subject to random packet dropouts which motivates us to consider a topology with random links. Given a network topology G we model the packet drops over unreliable links at time k by removing those links from the graph structure G(k). Simply speaking, given m packet drops over distinct links at time k, m links are removed from the network G(k). Then, to satisfy Assumption 2 for convergence, the remaining reliable network needs to hold uniform connectivity over every B time iterations for a B ∈ Z ≥0 .
Assuming common knowledge, both nodes i, j are aware of the delivery or loss of the messages or packets over the bi-directional link (i, j) [11], [13]. This is to (i) keep the networking balanced following Assumption 2 and (ii) satisfy the all-time feasibility in the proof of Theorem 1. For this purpose, in case the message from i to j is lost at iteration k, node i does not incorporate possibly received message g l (∂f j (x j (k)) from node j in its data-processing and updating state x i (k + 1) via (2). In other words, the mutual messages are either both dropped or both used in Algorithm ??. This consideration makes the probability of packet drop different from the probability of link removal in our calculations; for packet drop rate p d (or packet delivery rate 1−p d ) over the links from i to j or j to i, the equivalent probability of link removal in our analysis follows as implying the probability that either of the messages or both are lost.
Recall that the notion of uniform connectivity implies that the union network over every finite number of time-steps B is connected, i.e., G B = k+B k G(k) is connected. This is much more relaxed than the all-time connectivity requirement in many literature [1]- [8]. Intuitively, for example with B = 2, one can assume existence of B = 2 links between every two nodes i, j over the union network G B = k+B k G(k). Assume p l as the probability of link removal over connected graph G. Recalling the definition of the union graph, the link between i, j over G B (k) is lost if the link over both G(k) and G(k + 1) are removed as unreliable links. Therefore, the probability that any link is unreliable over the union network G B (k) is p 2 l for B = 2. One can extend this to uniform connectivity over any B iterations, i.e., the probability of link removal over This is the intuition behind the following theorem.
Theorem 2: Assume a connected network topology G with bond-percolation threshold p c . There exists B ∈ Z ≥0 such that G B = ∪ k+B k G(k) remains uniformly connected with probability 1 under any drop rate For p d < 1 − √ 1 − p c the network G(k) remains connected at iteration k with probability 1.
Proof: First, recall that for packet drop rate p d the probability of link removal is p l = 2p d − p 2 d from (9). For percolation threshold p c , find B * ∈ Z ≥0 as the minimum value satisfying p 1+B * l < p c . Note that, for p l < 1, function p 1+B * l is monotonically increasing on p l and decreasing on B * . Then, for any B ≥ B * , the probability that all the links are dropped over B * time-steps is less than the percolation threshold p c . This means that, for B ≥ B * and p c < p l < 1, the network G B remains uniformly connected with probability 1 under link removal probability p l = 2p d − p 2 d . This gives the admissible range of the packet drop rate 1 − √ 1 − p c < p d < 1. Corollary 1: Assume a dynamic sparse network which is not connected but uniformly-connected over B 0 > 0 time iterations. Given the bond-percolation threshold p c associated with G B0 and drop rate 1 d . Then, the union network G B under the same packet drop rate remains uniformly connected for any B ≥ B 0 B * .
Remark 2: In the case of heterogeneous and time-varying drop rates p ij (k) at different links (i, j), one can consider p d = max{p ij (k)} in Theorem 2 as a conservative solution to find minimum B * value.

V. SIMULATION
For the simulation we consider a network G of n = 20 nodes based on the ER model shown in Fig. 1(TopLeft). The linking probability between every two nodes is considered as p = 0.3. The ER theory states that for p > 1 n−1 , the network is connected [21], i.e., the connectivity transition occurs at We use MATLAB randperm function to randomly assign these rates over successive periods of every 200 steps. Sample network topologies (after removal of unreliable links) associated with the given probabilities p l are shown in Fig. 1. For the ER graph in Fig. 1(a) (with p l = 0), we have N = 5.6 (as average node degree) and, from Table I, we have p c = 0.177. Note that all considered p l values are over this threshold. From Theorem 2, we have respective values B * = [3,5,7,11,23], which implies uniform-connectivity over B ≥ 23 steps and, thus, any sufficiently small η satisfying Theorem 1 ensures convergence of the allocation algorithm, as described next.
To verify, we apply Algorithm ?? over this dynamic network to solve the following allocation problem: with random parameters a i , l i , c i , d i and adding penalty terms max{x i − M i , 0} 2 + max{m i − x i , 0} 2 to address the box constraints m i = 2, M i = 7. In the proposed dynamics (2), we consider two example strongly sign-preserving nonlinear functions g n (z) = z + z 3 and logarithmic quantizer g l (z) with ρ = 1 256 . We consider random symmetric link weights W ij in the range (0, 10] (non-stochastic). For the states bounded by the given box constraints, we have κ n = 1 and K n = 147 and u = 0.05. For the logarithmic quantization we have κ l = 1 − ρ 2 = 0.998 and K l = 1 + ρ 2 = 1.002. For the given network G B we have λ2 λ 2 n = 0.019. This gives η = 0.0025 as the bound on the step rate for convergence. This gives a sufficient bound to ensure convergence and, for faster decay rates, in this simulation we choose η = 0.05. Fig. 2 shows the evolution of the residual cost F and states x i in this example.

VI. DISCUSSIONS AND CONCLUDING REMARKS
This paper provides a robust approach for distributed resource allocation over networks with unreliable links that resemble communication channels subject to packet drops. By relaxing (i) the all-time network connectivity requirement to uniform connectivity and (ii) the dynamic weightstochastic condition to weight-balanced links, the proposed allocation algorithm is proved to converge under different packet drop rates, but over a longer time horizon. This approach finds application in mobile multi-agent networks and wireless communication networks with the inevitable high rate of packet losses. The proposed algorithm has many other advantages over the existing allocation solutions. For example, it can address nonlinearities on the agents' dynamics. This nonlinear model may resemble (i) physics of the system and inherent constraints in application, e.g., quantized information exchange [6], [9], clipping, or actuator saturation, and (ii) purposely designed dynamics to suppress impulsive noise or improve the convergence time [31]. The solution is not limited to the quadratic cost models, e.g., in CPU scheduling and economic dispatch, but it can handle general non-quadratic cost models, e.g., due to additive barrier functions [2] and penalty terms [26], [27]. As a future research direction, we are extending these results to address possible time-delays in data exchange over the links.
Another interesting preventive approach is survivable network design via Q-edge-connected graphs which is robust to link removal. This problem aims to design the network such that it remains connected after removing any subset of size (up to) Q links, or it preserves a prescribed routing criterion (up to) a certain disruption cost [34], [35]. By such a design, dropping up to Q packets over the network, a route/path of package delivery between every two nodes is guaranteed over time, ensuring convergence of our allocation algorithm. Recall that existing survivable design algorithms ensure connectivity at every time-instant and extension to more relaxed uniform connectivity over time (as described in Assumption 2) is another promising future research direction.