s2fieldsofstudy:
Computer Science
Economics

text:
Comparative Analysis of Financial Network Topology for the Russian, Chinese and US Stock Markets

This paper studies the properties of the Russian stock market by employing the data-driven science and network approaches. The theory of complex networks allows us to build and examine topological network structures of the market with the further identification of relationships between stocks and the analysis of hidden information and market dynamics. In this paper we will present an analysis of structural and topological properties of the Russian stock market using market graph, hierarchical tree, minimum spanning tree approaches. We compare topological properties of the networks constructed for the US and China stock markets with the properties of correspondingnetworksconstructedfortheRussianstockmarketusingadatasetspanningovereightyears. Key–Words: Network analysis, Market graph, Degree distribution, Hierarchical tree, Minimum spanning tree Received: November 7, 2019. Revised: March 9, 2020. Accepted: March 25, 2020. Published: March 31, 2020.

Introduction
To represent the interaction framework between agents in economic or financial systems, one can use network models. The definition of a financial network was given in [1], in which the nodes are various organizations, including banks, firms and investors, which are linked through a network of financial interdependencies, for example, cross-ownership of financial assets, inter-organizational debts / liabilities, social relations between board members and so on. A classic example of a banking (financial) network is built in [2], where it is illustrated how systemic risk can arise as a result of financial interdependencies be-tween banks. Most of real-world networks are constantly changing over time. The links between economic agents are usually activated or deactivated in dependence of the irregular interactions between the elements of the economic system. These activation templates may be sporadic and represent a huge role in the dynamics of network processes. One of examples of this situation is information spreading in social networks. The topology of network (heterogeneity in the connections, abundance of loops) plays an important role in the economic system dynamics behavior.
Financial markets and corresponding networks may have a huge amount of assets (nodes) and inter-connections (edges). It is useful to simplify initial networks by filtering the noises. Some of the approaches employed in filtering of financial networks are based on hierarchical tree (HT) [3,4], planar maximally filtered graphs (PMFG) [5], asset graph (AG) [6], and partial correlation network [7,8,9]. Another successful method of financial network simplification was developed by Mantegna [3] and applied in [10]. The approach is based on the minimum spanning tree method for examining of hierarchical structures of financial networks.
In recent years there has been an increased interest in the application and development of an approach based on market graphs. These papers include empirical studies based on real market data and explore various structural properties and attributes of market graphs, such as maximum clicks, maximum independent sets, degree distribution [11,12,13,14], clustering, Pearson correlations [15], market graphs of the US market [16], the complexity of the market graph [17]. The concept of a market graph was first considered in [18], in which a market network is defined as a full weighted graph, where nodes represent assets and arc weights reflect similarities between asset behavior (which can be measured, for example, using correlation). In article [18], the edge between two vertices is inserted into the market graph, if the corresponding value of the correlation coefficient is above a given threshold. The articles [19,20,21,12,22] studied the distinctive features of individual financial markets. Market graphs with similarity measures other than correlation were studied in [19,23,24,25,8,26].
In this paper, we use some methods of complex network theory to study Russian financial market network. Complex network analysis and SNA methods allow us to investigate the structure of Russian financial market network that are based on the correlation relationship between the components of the system (assets).
It turned out that the network analysis is capable to examine the importance of companies using different centrality measures [45,45] and to analyze the systemic risks and financial market stability based on the topological properties of networks [46,47,48]. Systemic risk contagion has been analyzed in [49,50] to find major influencers based on banking behaviors in the global environment. The paper [51] studies European markets and shows that EU markets are exposed to systemic risks. It should be noted that systemic risk analysis based on network approach may provide useful implications for market regulations.
The paper [29] employs various methods of network analysis to investigate financial networks of both in China and the United States and to study how the two markets behave differently. In this paper, we will present a similar analysis of Russian stock market. One of the main research questions of this paper is to find how the US and China stock markets differ in the structures and topological properties from Russian stock market using • market graph approach proposed in [18,11], • hierarchical tree (HT) method [3], • minimum spanning tree (MST) approach [3].
In this paper we compare topological properties of the networks constructed for the US and China stock markets in [29] with the properties of corresponding networks constructed for the Russian stock market using a dataset spanning over nine years.
The paper is organized as follows. First, Section 2 describes the data and method used to construct networks. Section 3 presents the market network properties using market graph approach proposed in [18,11]. Section 4 examines the hierarchical structures of the Russian stock network using HT, and MST approaches and compares them with findings of the paper [29]. Finally, conclusions are presented in Section 5.

Data Description
Russia has a developing economy which is the sixth among countries in the world in terms of GDP (PPP). The RTS index and the Moscow Stock Exchange index are the main indicators of the Russian stock market. The difference between them is that the RTS is calculated in dollars, and the MICEX in rubles. Historical data of daily prices for each company were taken from open sources on Yahoo Finance, for the period from 10/01/2012 to 09/04/2019 (1664 trading days). Similarly, data on the MICEX index were taken from the Moscow Stock Exchange for the same period. A total of 194 companies are considered, among which 32 companies are part of the MICEX index.
The descriptive statistics for log returns of the IMOEX index over the period between 10/01/2012 and 09/04/2019 are presented in Table 2.2. The result of Garch(1,1) modelling for log returns of the IMOEX index over the period between 10/01/2012 and 09/04/2019 is presented in Fig. 1.
Some companies (CBOM, DSKY, FIVE, LNTA, MOEX, POLY, RNFT, RUAL, SFIN) that are part of the MICEX index were not included in our dataset. This is due to the fact that they have a very short trading history (only 2-3 years), which is insufficient to build a relevant model. We also did not include companies that had more than 250 missed trading days in a row in our dataset. Thus, only 194 issuers remained in our dataset from the initial list of 278 issuers traded on the Moscow stock exchange. The YNDX company has started to be traded on the Moscow stock exchange from 02/06/2014, but it has been traded on the NASDAQ from 02/05/2011, from which the missing data were taken. After cleaning and the selection process, we include 194 stocks traded on IMOEX in our dataset, as shown in Table 1. We place these 194 stocks among all ten industry sectors.

Market Network Construction
To calculate the correlation ρ ij for a pair of stocks, it is necessary to use time series of prices (Adj Close) P i (t) for each company s i at the time point t. To smooth the oscillations we use the log returns Y i (t) of the company s i in the time where ∆t = 1 for daily prices. Then we find the Pearson correlation coefficient for each pair of companies s i and s j as follows where · denotes the average value.
Based on the values of the Pearson correlation coefficient for each pair of companies we can construct the N × N = 194 × 194 distance matrix using the equation

Basic network properties
In this section, we examine the topological properties of the IMOEX stock network using the dynamic approach applied in [29]. We take the sliding window size L = 250 to meet the requirement of L/N > 1, where the number of stocks is N = 194. In total, there are 1413 windows for IMOEX. We calculated the log returns for all IMOEX stocks by using Eq.
(1) and then we found the Pearson correlation coefficient matrix for each window over the period between 10/01/2012 and 09/04/2019 using Eq. (2). We compare the results obtained for the Russian stock market with the results obtained in the paper [29] for the Chinese and US stock markets. The dataset S&P468 has been described in [29] Table 2: Basic information of IMOEX datasets including the number of stocks, the number of trading days, the average log returns r , and the standard deviation σ r of the log returns is presented dataset CSI163 has been described in [29] and includes of 163 component stocks of CSI300. Thus, for IMOEX we obtain a network with 194 vertices. Using Eq. (3) we calculate the distance matrices for each window. It should be noted that all negative correlation values are transformed into positive distances. All vertices in the network are fixed and correspond to stocks. However, the edges vary in each sliding window as the distances change. Then we average all the distance matrices over all sliding windows of the study periods. The statistical properties of the IMOEX stock network are examined in this section.
The maximum number of edges in average distance matrix is equal to N (N − 1)/2 = 18, 721 and it is a huge number. We can simplify the network by filtering edges with small weights. In this section the market graph is constructed based on threshold approach: for a fixed threshold value of θ we delete those edges whose weight are bigger than θ and preserve the remaining edges. In other words, we filter the network edges as follows: Table 3 presents the basic properties of the IMOEX network for different values of θ taken from 0.5 to 1.5. It can be seen that the number of edges is zero if the threshold value is θ = 0.8, while the graph becomes complete at θ = 1.5. The figure 2 shows edge densities for different thresholds from 0.5 to 1.5. In the interval from 0.5 to 1.1, network densities are close to 0, which means that all edges are filtered. In the range from 1.43 to 1.5, the densities are close to 1, i.e. all edges are present at the network. For charts based on Chinese and American stock markets [29] we see that the two curves have a similar shape with a slopes that lie between 0.6 and 1. It should be noted that the transformation interval for the IMOEX network is from 1.1 to 1.43, which is atypical for stock networks. As a rule, stock networks show the transformation interval from 0.6 to 1. The study shows that Russian stock network has a different form of trans-forming interval than the USA and China . This is due to much lower correlations between assets traded in the Russian stock market, which could be explained by e.g. low liquidity of most stocks in IMOEX. Moreover, it can be easily verified that even highest liquidity shares provide low correlations. Network density for a graph constructed using these stocks maintains similar dependency from θ as for the whole IMOEX network.

Degree distribution
We investigate the degree distributions for the IMOEX networks filtered with different values of θ. The degree distributions are noisy for both small values of θ and large ones. For a sufficiently small interval θ = (1.15; 1.2), the degree distributions follow closely to the power law distribution. After processing the data, we build the degree distribution with the typical power law shape in Fig. 3. The log-binning procedure was used to build the regression and calculate the power law exponent. The value of the power law exponent is γ = 1.14 with R 2 = 0.85. Since the use of linear approximation (linear-binning) significantly underestimates the exponent γ. A negative tilt angle indicates that the network is scale-free, in which a small part of  are presented for different θ from 0.5 to 1.5 in a step of 0.1, the standard deviation σ the vertices has greater degrees, and most of the vertices have smaller degrees. Fig. 3 shows that the degree distributions of IMOEX (with θ = 1.2), CSI163 (with θ = 0.68) and S&P468 (with θ = 0.75) networks follow the power law in the form of p k ∼ k −γ with γ = 1.14, γ = 0.9935 and γ = 1.2323 respectively. It means that the networks are scale free in which a small portion of vertices have larger degrees, while a large portion of vertices have small degrees. At the first sight, the IMOEX network has the same scale-free behaviour as the CSI163 and S&P468 networks. However, if we take the same level of θ as it was taken for the CSI163 and S&P468 networks (i.e. θ = 0.65 ± 0.1) then the IMOEX network will be very sparse and almost all vertices will be isolated.

Average clustering coefficient
The average clustering coefficient (ACC) is the average of all clustering coefficients over all vertices. The clustering coefficient means transitivity for a single vertex, while the average clustering coefficient is an indicator of the transitivity and density for the entire network. Fig. 4 presents the average clustering coefficient for the IMOEX network and its comparison with a random network. Compared to the random network of the same size, the ACC for the IMOEX network is significantly larger than that of the random network. It shows that the IMOEX network is well connected with better transitivity. For comparison, Fig. 4 presents the average clustering coefficients for China and the USA networks. Fig. 4 shows that the IMOEX network does not behave in the similar way: ACC values are always lesser than for CSI163 and S&P468 networks, i.e. the transitivity and density for the IMOEX network are much lower than for China and the USA networks.

Average path length
Fig . 5 shows the average path length for the IMOEX network compared to a random network of the same size. It can be seen that the IMOEX network is significantly different from the random network. For comparison, Fig. 5 also shows the average clustering coefficients for China and the USA networks. Fig.  5 indicates that the IMOEX network behaves differently: the average path length has a peak at θ = 1.1 while the CSI163 and S&P468 networks have peaks at θ = 0.6 and 0.65 respectively. Moreover, the average path length for the IMOEX network is much bigger than for the CSI163 and S&P468 networks for θ > 1. Thus, the China and the USA networks are more con-

Betweenness centrality
Betweenness of a vertex is defined as the number of shortest paths passing through it. Thus, this indicator demonstrates the contribution of each vertex to global connectivity. Averaging over individual betweenness coefficients over all vertices, we can compare the betweenness values for any two vertices. Fig. 6 shows that the peak of the average betweenness for the network is achieved at the threshold value θ = 1.35. It can be seen (Fig. 5)   We visualize the stock network of IMOEX with different values of θ = 1.1, 1.2 in Fig. 7. The figures show that the network can be greatly simplified by using small values of θ. Table 3 indicates that the edge density for the IMOEX network increases dramatically from θ = 1.1 to θ = 1.4. To examine how industry sectors are linked in the IMOEX network, we list the properties of the network with θ = 1.4 in Table 3.5. The table shows that the Russian industry sectors has approximately the same average degree as well as the average clustering coefficient. However, the values of average betweenness coefficient differ for different sectors assuredly. It should be noted that the same results were obtained for the CSI163 and S&P468 networks in [29].   Table 4: In this table, we list the sector name and the numbers of stocks, the average degree d , the average clustering coefficient C , and the average betweenness coefficient B for each industry sector of these 194 stocks. The values are calculated from the IMOEX network with θ = 1.4. Table 3.5 presents the top 10 stocks of IMOEX network with the biggest values of degree d i and betweenness b i arranged in descending order in the top part and the bottom part, respectively. The stock code, company name, industry name, and values of degree d i , clustering coefficient c i and betweenness b i are presented. SBERBANK, which is a leading financial company in Russia, has the largest degree of 136 and JSC "OGK-2" is Russia's largest heat generating company, has the largest betweenness coefficient of 1551.

Hierarchical Structure of Stock Networks
Graphs constructed by certain rules may reveal its properties. In this section we create and study hierarchical trees and minimal spanning tree to extract valuable fragments from the graph based on the average distance matrix, which we construct by averaging sliding windows over the study period.

Hierarchical Tree
Clustering algorithms can be used on average distance matrix to group stocks which prices correlate and to   malised distance definition (Fig. 11). Dendogram function from Scipy package is used to draw them on a plot from Matplotlib package. In the trees, each leaf node indicates the stock of a particular company with its ticket as a label. The height of where two stocks connect shows how they correlate. The lower the branches of similar clusters or stocks merge the more positive correlation there is between them. We set color threshold of 0.7×max(d * ) to highlight clusters below this value. The companies in same colored branches are very likely to be from the same economy

Minimum Spanning Tree
IMOEX network, based on average distance matrix, can be simplified by extracting edges with minimal distance while maintaining graph connectivity. Connected graph with N vertices can be converted into the graph with N − 1 edges without loops with minimum total length of edges or minimum spanning tree (MST). This graph is called minimum spanning tree. It brings huge advantages to the study of networks of stocks by reducing noises and simplifying computations. We use Python Scipy minimum spanning tree algorithm on an average distance matrix. Networkx python package is used to get its graphic representation on a Matplotlib plot as shown in Fig. 12.
After the edge filtering process, some nodes are still well-connected, and they may be considered the key stocks in the market. The most connected stocks are from utilities, financials and energy sectors. These results let us extract the most influential companies from the given data. It should be noted that Markowitz optimal portfolio consists of stocks corresponding to leaf nodes of MST. Therefore, we can use MST approach for portfolio optimisation.
Top 10 stocks are shown in the table 4.2. As expected, financials, energy and utilities sectors are dominant in the Russian stock market which is different from the US and Chinese stock markets. Further studies of minimum spanning tree over time may reveal important changes in the core stocks and sectors of the market. Moreover, the result of the algorithm application differs between study periods and this is helpful when studying network topological dynamics.

Hierarchical tree for a pooled market graph (IMOEX and S&P468)
The obtained results indicate a significant difference in the density of correlation between stock returns and the characteristics of market graphs constructed for the financial markets of the US, China and the Russian Federation. Such differences can be explained by the relatively less developed financial market of the Russian Federation as a whole, and industry affiliation of companies, greater uncertainty and a high share of risks specific to each of the companies, etc. It is known that a significant share in the total capitalization of the financial market of the Russian Federation falls   on the shares of companies focused on the export of raw materials. It is interesting to check whether industry or country is more important for the correlation dependence. To do this, we apply the hierarchical tree method for a pooled market graph (IMOEX and S&P468). If industry factors would be more important, then Russian stocks should be placed to clusters according to their industry affiliation. If cross-country differences are more important, then clusters should be placed to clusters according to their country affili-ation.
The results of applying the hierarchical classification method are presented in the Fig. 13. All Russian companies fall into one cluster (in the figure it is marked in blue at the left side of the dendrogram). S&P500 companies are located in other clusters, which are well explained by their industry affiliation.

Conclusion
In this paper, using the daily close prices we calculated the correlation matrices for 194 assets trading on the Russian stock exchange during the study period from 2012 to 2019 for all sliding windows. Based on these correlation matrices we constructed the market networks with stocks as the nodes and correlation relationships as the edges. Our goal was to study the networks properties using the complex network theory from data science perspective. We employed such network filtering methods as threshold, hierarchical tree, minimum spanning tree for simplification of the networks. Then we examine the network properties for the filtered networks. We compare topological properties of the networks constructed for the US and China stock markets with the properties of corresponding networks constructed for the Russian stock market using a dataset spanning over eight years. It turns out that the IMOEX network share similar topological properties in some case, but it also differ from the CSI163 and S&P468 networks in many aspects. The study shows that the topological properties of the Russian stock network has a different from the USA and China dependence on θ, which is the consequence of a much lower correlations between assets traded in the Russian stock market. One the one hand, it may be explained by low liquidity of most stocks traded in IMOEX. On the other hand, even shares with high liquidity have low correlations.


extfieldsofstudy:
Computer Science
Engineering
Mathematics

text:
SPACE: Unsupervised Object-Oriented Scene Representation via Spatial Attention and Decomposition

The ability to decompose complex multi-object scenes into meaningful abstractions like objects is fundamental to achieve higher-level cognition. Previous approaches for unsupervised object-oriented scene representation learning are either based on spatial-attention or scene-mixture approaches and limited in scalability which is a main obstacle towards modeling real-world scenes. In this paper, we propose a generative latent variable model, called SPACE, that provides a unified probabilistic modeling framework that combines the best of spatial-attention and scene-mixture approaches. SPACE can explicitly provide factorized object representations for foreground objects while also decomposing background segments of complex morphology. Previous models are good at either of these, but not both. SPACE also resolves the scalability problems of previous methods by incorporating parallel spatial-attention and thus is applicable to scenes with a large number of objects without performance degradations. We show through experiments on Atari and 3D-Rooms that SPACE achieves the above properties consistently in comparison to SPAIR, IODINE, and GENESIS. Results of our experiments can be found on our project website: https://sites.google.com/view/space-project-page

INTRODUCTION
One of the unsolved key challenges in machine learning is unsupervised learning of structured representation for a visual scene containing many objects with occlusion, partial observability, and complex background. When properly decomposed into meaningful abstract entities such as objects and spaces, this structured representation brings many advantages of abstract (symbolic) representation to areas where contemporary deep learning approaches with a global continuous vector representation of a scene have not been successful. For example, a structured representation may improve sample efficiency for downstream tasks such as a deep reinforcement learning agent (Mnih et al., 2013). It may also enable visual variable binding (Sun, 1992) for reasoning and causal inference over the relationships between the objects and agents in a scene. Structured representations also provide composability and transferability for better generalization.
Recent approaches to this problem of unsupervised object-oriented scene representation can be categorized into two types of models: scene-mixture models and spatial-attention models. In scenemixture models (Greff et al., 2017;Engelcke et al., 2019), a visual scene is explained by a mixture of a finite number of component images. This type of representation provides flexible segmentation maps that can handle objects and background segments of complex morphology. However, since each component corresponds to a full-scale image, important physical features of objects like position and scale are only implicitly encoded in the scale of a full image and further disentanglement is required to extract these useful features. Also, since it does not explicitly reflect useful inductive biases like the locality of an object in the Gestalt principles (Koffka, 2013), the resulting component representation is not necessarily a representation of a local area. Moreover, to obtain a complete scene, a component needs to refer to other components, and thus inference is inherently performed sequentially, resulting in limitations in scaling to scenes with many objects.
In contrast, spatial-attention models (Eslami et al., 2016;Crawford & Pineau, 2019) can explicitly obtain the fully disentangled geometric representation of objects such as position and scale. Such features are grounded on the semantics of physics and should be useful in many ways (e.g., sample efficiency, interpretability, geometric reasoning and inference, transferability). However, these models cannot represent complex objects and background segments that have too flexible morphology to be captured by spatial attention (i.e. based on rectangular bounding boxes). Similar to scene-mixture models, previous models in this class show scalability issues as objects are processed sequentially.
In this paper, we propose a method, called Spatially Parallel Attention and Component Extraction (SPACE), that combines the best of both approaches. SPACE learns to process foreground objects, which can be captured efficiently by bounding boxes, by using parallel spatial-attention while decomposing the remaining area that includes both morphologically complex objects and background segments by using component mixtures. Thus, SPACE provides an object-wise disentangled representation of foreground objects along with explicit properties like position and scale per object while also providing decomposed representations of complex background components. Furthermore, by fully parallelizing the foreground object processing, we resolve the scalability issue of existing spatial attention methods. In experiments on 3D-room scenes and Atari game scenes, we quantitatively and qualitatively compare the representation of SPACE to other models and show that SPACE combines the benefits of both approaches in addition to significant speed-ups due to the parallel foreground processing.
The contributions of the paper are as follows. First, we introduce a model that unifies the benefits of spatial-attention and scene-mixture approaches in a principled framework of probabilistic latent variable modeling. Second, we introduce a spatially parallel multi-object processing module and demonstrate that it can significantly mitigate the scalability problems of previous methods. Lastly, we provide an extensive comparison with previous models where we illustrate the capabilities and limitations of each method.

THE PROPOSED MODEL: SPACE
In this section, we describe our proposed model, Spatially Parallel Attention and Component Extraction (SPACE). The main idea of SPACE, presented in Figure 1, is to propose a unified probabilistic generative model that combines the benefits of the spatial-attention and scene-mixture models.

GENERATIVE PROCESS
SPACE assumes that a scene x is decomposed into two independent latents: foreground z fg and background z bg . The foreground is further decomposed into a set of independent foreground objects z fg = {z fg i } and the background is also decomposed further into a sequence of background segments z bg = z bg 1:K . While our choice of modeling the foreground and background independently worked well empirically, for better generation, it may also be possible to condition one on the other. The image distributions of the foreground objects and the background components are combined together with a pixel-wise mixture model to produce the complete image distribution: Here, the foreground mixing probability α is computed as α = f α (z fg ). This way, the foreground is given precedence in assigning its own mixing weight and the remaining is apportioned to the background. The mixing weight assigned to the background is further sub-divided among the K background components. These weights are computed as π k = f π k (z bg 1:k ) and k π k = 1. With these notations, the complete generative model can be described as follows.
p(x) = p(x|z fg , z bg )p(z bg )p(z fg )dz fg dz bg (2) Figure 1: Illustration of the SPACE model. SPACE consists of a foreground module and a background module.
In the foreground module, the input image is divided into a grid of H × W cells (4 × 4 in the figure). An image encoder is used to compute the z where , z depth , and z pres for each cell in parallel. z where is used to identify proposal bounding boxes and a spatial transformer is used to attend to each bounding box in parallel, computing a z what encoding for each cell. The model selects patches using the bounding boxes and reconstructs them using a VAE from all the foreground latents z fg . The background module segments the scene into K components (4 in the figure) using a pixel-wise mixture model. Each component consists of a set of latents z bg = (z m , z c ) where z m models the mixing probability of the component and z c models the RGB distribution of the component. The components are combined to reconstruct the background using a VAE. The reconstructed background and foreground are then combined using a pixel-wise mixture model to generate the full reconstructed image.
We now describe the foreground and background models in more detail.
Foreground. SPACE implements z fg as a structured latent. In this structure, an image is treated as if it were divided into H × W cells and each cell is tasked with modeling at most one (nearby) object in the scene. This type of structuring has been used in (Redmon et al., 2016;Santoro et al., 2017;Crawford & Pineau, 2019). Similar to SPAIR, in order to model an object, each cell i is associated with a set of latents (z pres In this notation, z pres is a binary random variable denoting if the cell models any object or not, z where denotes the size of the object and its location relative to the cell, z depth denotes the depth of the object to resolve occlusions and z what models the object appearance and its mask. These latents may then be used to compute the foreground image component p(x|z fg ) which is modeled as a Gaussian distribution N (µ fg , σ 2 fg ). In practice, we treat σ 2 fg as a hyperparameter and decode only the mean image µ fg . In this process, SPACE reconstructs the objects associated to each cell having z pres i = 1. For each such cell, the model uses the z what i to decode the object glimpse and its mask and the glimpse is then positioned on a full-resolution canvas using z where i via the Spatial Transformer (Jaderberg et al., 2015). Using the object masks and z depth i , all the foreground objects are combined into a single foreground mean-image µ fg and the foreground mask α (See Appendix D for more details).
SPACE imposes a prior distribution on these latents as follows: Here, only z pres i is modeled using a Bernoulli distribution while the remaining are modeled as Gaussian.
Background. To model the background, SPACE implements z bg k , similar to GENESIS, as (z m k , z c k ) where z m k models the mixing probabilities π k of the components and z c k models the RGB distribution p(x|z bg k ) of the k th background component as a Gaussian N (µ bg i , σ 2 bg ). The following prior is imposed upon these latents.

INFERENCE AND TRAINING
Since we cannot analytically evaluate the integrals in equation 2 due to the continuous latents z fg and z bg 1:K , we train the model using a variational approximation. The true posterior on these variables is approximated as follows.
This is used to derive the following ELBO to train the model using the reparameterization trick and SGD (Kingma & Welling, 2013).
See Appendix B for the detailed decomposition of the ELBO and the related details.
Parallel Inference of Cell Latents. SPACE uses mean-field approximation when inferring the cell latents, so z fg for each cell does not depend on other cells.
As shown in Figure 1, this allows each cell to act as an independent object detector, spatially attending to its own local region in parallel. This is in contrast to inference in SPAIR, where each cell's latents auto-regressively depend on some or all of the previously traversed cells in a row-major order i.e., q(z fg |x) = HW i=1 q(z fg i |z fg <i , x). However, this method becomes prohibitively expensive in practice as the number of objects increases. While Crawford & Pineau (2019) claim that these lateral connections are crucial for performance since they model dependencies between objects and thus prevent duplicate detections, we challenge this assertion by observing that 1) due to the bottom-up encoding conditioning on the input image, each cell should have information about its nearby area without explicitly communicating with other cells, and 2) in (physical) spatial space, two objects cannot exist at the same position. Thus, the relation and interference between objects should not be severe and the mean-field approximation is a good choice in our model. In our experiments, we verify empirically that this is indeed the case and observe that SPACE shows comparable detection performance to SPAIR while having significant gains in training speeds and efficiently scaling to scenes with many objects.
Preventing Box-Splitting. If the prior for the bounding box size is set to be too small, then the model could split a large object by multiple bounding boxes and when the size prior is too large, the model may not capture small objects in the scene, resulting in a trade-off between the prior values of the bounding box size. To alleviate this problem, we found it sometimes helpful to introduce an auxiliary loss which we call the boundary loss. In the boundary loss, we construct a boundary of thickness b pixels along the borders of each glimpse. Then, we restrict an object to be inside this boundary and penalize the model if an object's mask overlaps with the boundary area. Thus, the model is penalized if it tries to split a large object by multiple smaller bounding boxes. A detailed implementation of the boundary loss is mentioned in Appendix C.

RELATED WORKS
Our proposed model is inspired by several recent works in unsupervised object-oriented scene decomposition. The Attend-Infer-Repeat (AIR) (Eslami et al., 2016) framework uses a recurrent neural network to attend to different objects in a scene and each object is sequentially processed one at a time. An object-oriented latent representation is prescribed that consists of 'what', 'where', and 'presence' variables. The 'what' variable stores the appearance information of the object, the 'where' variable represents the location of the object in the image, and the 'presence' variable controls how many steps the recurrent network runs and acts as an interruption variable when the model decides that all objects have been processed.
Since the number of steps AIR runs scales with the number of objects it attends to, it does not scale well to images with many objects. Spatially Invariant Attend, Infer, Repeat (SPAIR) (Crawford & Pineau, 2019) attempts to address this issue by replacing the recurrent network with a convolutional network. Similar to YOLO (Redmon et al., 2016), the locations of objects are specified relative to local grid cells rather than the entire image, which allow for spatially invariant computations. In the encoder network, a convolutional neural network is first used to map the image to a feature volume with dimensions equal to a pre-specified grid size. Then, each cell of the grid is processed sequentially to produce objects. This is done sequentially because the processing of each cell takes as input feature vectors and sampled objects of nearby cells that have already been processed. SPAIR therefore scales with the pre-defined grid size which also represents the maximum number of objects that can be detected. Our model uses an approach similar to SPAIR to detect foreground objects, but importantly we make the foreground object processing fully parallel to scale to large number of objects without performance degradation. Works based on Neural Expectation Maximization (Van Steenkiste et al., 2018;Greff et al., 2017) do achieve unsupervised object detection but do not explicitly model the presence, appearance, and location of objects. These methods also suffer from the problem of scaling to images with a large number of objects. In a related line of research, AIR has also recently been extended to track objects in a sequence of images (Kosiorek et al., 2018;Crawford & Pineau, 2020;Jiang et al., 2020).
For unsupervised scene-mixture models, several recent models have shown promising results. MONet (Burgess et al., 2019) leverages a deterministic recurrent attention network that outputs pixel-wise masks for the scene components. A variational autoencoder (VAE) (Kingma & Welling, 2013) is then used to model each component. IODINE (Greff et al., 2019) approaches the problem from a spatial mixture model perspective and uses amortized iterative refinement of latent object representations within the variational framework. GENESIS (Engelcke et al., 2019) also uses a spatial mixture model which is encoded by component-wise latent variables. Relationships between these components are captured with an autoregressive prior, allowing complete images to be modeled by a collection of components.

EVALUATION
We evaluate our model on two datasets: 1) an Atari (Bellemare et al., 2013) dataset that consists of random images from a pretrained agent playing the games, and 2) a generated 3D-room dataset that consists of images of a walled enclosure with a random number of objects on the floor. In order to test the scalability of our model, we use both a small 3D-room dataset that has 4-8 objects and a large 3D-room dataset that has 18-24 objects. Each image is taken from a random camera angle and the colors of the objects, walls, floor, and sky are also chosen at random. Additional details of the datasets can be found in the Appendix E.
Baselines. We compare our model against two scene-mixture models (IODINE and GENESIS) and one spatial-attention model (SPAIR). Since SPAIR does not have an explicit background component, we add an additional VAE for processing the background. Additionally, we test against two implementations of SPAIR: one where we train on the entire image using a 16 × 16 grid and another where we train on random 32 × 32 pixel patches using a 4 × 4 grid. We denote the former model as SPAIR and the latter as SPAIR-P. SPAIR-P is consistent with the SPAIR's alternative training regime on Space Invaders demonstrated in Crawford & Pineau (2019) to address the slow training of SPAIR on the full grid size because of its sequential inference. Lastly, for performance reasons, unlike the original SPAIR implementation, we use parallel processing for rendering the objects from their respective latents onto the canvas 1 for both SPAIR and SPAIR-P. Thus, because of these improvements, our SPAIR implementation can be seen as a stronger baseline than the original SPAIR. More details of the baselines are given in Appendix D.

QUALITATIVE COMPARISON OF INFERRED REPRESENTATIONS
In this section, we provide a qualitative analysis of the generated representations of the different models. For each model, we performed a hyperparameter search and present the results for the best settings of hyperparameters for each environment. Figure 2 shows sample scene decompositions of our baselines from the 3D-Room dataset and Figure 3 shows the results on Atari. Note that SPAIR does not use component masks and IODINE and GENESIS do not separate foreground from background, hence the corresponding cells are left empty. Additionally, we only show a few representative components for IODINE and GENESIS since we ran those experiments with larger K than can be displayed. More qualitative results of SPACE can be found in Appendix A.

IODINE & GENESIS.
In the 3D-Room environment, IODINE is able to segment the objects and the background into separate components. However, it occasionally does not properly decompose objects (in the Large 3D-room results, the orange sphere on the right is not reconstructed) and may generate blurry objects. GENESIS is able to segment the background walls, floor, and sky into multiple components. It is able to capture blurry foreground objects in the Small 3D-Room, but is not able to cleanly capture foreground objects with the larger number of objects in the Large 3D-Room. In Atari, both IODINE and GENESIS fail to capture the foreground properly or try to encode all objects in a single component. We believe this is because the objects in Atari games are smaller, less regular and lack the obvious latent factors like color and shape as in the 3D dataset, and thus detection-based approaches are more appropriate in this case.
SPAIR & SPAIR-P. SPAIR is able to detect tight bounding boxes in both 3D-Rooms and the two Atari games. SPAIR-P, however, often fails to detect the foreground objects in proper bounding boxes, frequently uses multiple bounding boxes for one object and redundantly detects parts of the background as foreground objects. This is a limitation of the patch training as the receptive field of  each patch is limited to a 32 × 32 glimpse, prohibiting it from detecting larger objects and making it difficult to distinguish the background from foreground. These two properties are illustrated well in Space Invaders, where it is able to detect the small aliens, but it detects the long piece of background ground on the bottom of the image as foreground objects.

SPACE.
In both 3D-Room, SPACE is able to accurately detect almost all objects despite the large variations in object positions, colors, and shapes, while producing a clean segmentation of the background walls, ground, and sky. This is in contrast to the SPAIR model, while being able to provide similar foreground detection quality, encodes the whole background into a single component, which makes the representation less disentangled. Notably, in River Raid where the background is constantly changing, SPACE is able to perfectly segment the blue river while accurately detecting all foreground objects, while SPAIR often cannot properly separate foreground and background.
Joint Training. Figure 4 shows the results of training SPACE jointly across 10 Atari games. We see that even in this setting, SPACE is able to correctly detect foreground objects and cleanly segment the background.
Foreground vs Background. Typically, foreground is the dynamic local part of the scene that we are interested in, and background is the relatively static and global part. This definition, though intuitive, is ambiguous. Some objects, such as the red shields in Space Invaders and the key in Montezuma's Revenge ( Figure 6) are static but important to detect as foreground objects. We found that SPACE tends to detect these as foreground objects while SPAIR considers it background. Similar behavior is observed in Atlantis (Figure 4), where SPACE tends to detect some foreground objects from the middle base that is above the water. One reason for this behavior is because we limit the capacity of the background module by using a spatial broadcast decoder (Watters et al., 2019) which is much weaker when compared to other decoders like sub-pixel convolutional nets (Shi et al. (2016)). This would favor modeling static objects as foreground rather than background.

QUANTITATIVE COMPARISON
In this section we compare SPACE with the baselines in several quantitative metrics 2 . We first note that each of the baseline models has a different decomposition capacity (C), which we define as the capability of the model to decompose the scene into its semantic constituents such as the foreground objects and the background segmented components. For SPACE, the decomposition capacity is equal to the number of grid cells H × W (which is the maximum number of foreground objects that can be detected) plus the number of background components K. For SPAIR, the decomposition capacity is equal to the number of grid cells H × W plus 1 for background. For IODINE and GENESIS, it is equal to the number of components K.
For each experiment, we compare the metrics for each model with similar decomposition capacities. This way, each model can decompose the image into the same number of components. For a setting in SPACE with a grid size of H × W with K SPACE components, the equivalent settings in IODINE and GENESIS would be with C = (H × W ) + K SPACE . The equivalent setting in SPAIR would be a grid size of H × W .

Gradient
Step Latency. The leftmost chart of Figure 5 shows the time taken to complete one gradient step (forward and backward propagation) for different decomposition capacities for each of the models. We see that SPAIR's latency grows with the number of cells because of the sequential  nature of its latent inference step. Similarly GENESIS and IODINE's latency grows with the number of components K because each component is processed sequentially in both the models. IODINE is the slowest overall with its computationally expensive iterative inference procedure. Furthermore, both IODINE and GENESIS require storing data for each of the K components, so we were unable to run our experiments on 256 components or greater before running out of memory on our 22GB GPU. On the other hand, SPACE employs parallel processing for the foreground which makes it scalable to large grid sizes, allowing it to detect a large number of foreground objects without any significant performance degradation. Although this data was collected for gradient step latency, this comparison implies a similar relationship exists with inference time which is a main component in the gradient step.
Time for Convergence. The remaining three charts in Figure 5 show the amount of time each model takes to converge in different experimental settings. We use the pixel-wise mean squared error (MSE) as a measurement of how close a model is to convergence. In all settings, SPAIR and SPACE converge much faster than IODINE and GENESIS. In the 4 × 4 and 8 × 8 setting, SPAIR and SPACE converge equally fast. But as we scale up to 16 × 16, SPAIR becomes much slower than SPACE .
Average Precision and Error Rate. In order to assess the quality of our bounding box predictions, we measure the Average Precision and Object Count Error Rate of our predictions. Our results are shown in Table 1. We only report these metrics for 3D-Room since we have access to the ground truth bounding boxes for each of the objects in the scene. Both models have very similar average precision and error rate. Despite being parallel in its inference, SPACE has a comparable count error rate to that of SPAIR.
From our experiments, we can assert that SPACE can produce similar quality bounding boxes as SPAIR while 1) having orders of magnitude faster inference and gradient step time, 2) scaling to a large number of objects without significant performance degradation, and 3) providing complex background segmentation.

CONCLUSION
We propose SPACE, a unified probabilistic model that combines the benefits of the object representation models based on spatial attention and the scene decomposition models based on component mixture. SPACE can explicitly provide factorized object representation per foreground object while also decomposing complex background segments. SPACE also achieves a significant speed-up and thus makes the model applicable to scenes with a much larger number of objects without performance degradation. Besides, the detected objects in SPACE are also more intuitive than other methods. We show the above properties of SPACE on Atari and 3D-Rooms. Interesting future directions are to replace the sequential processing of background by a parallel one and to improve the model for natural images. Our next plan is to apply SPACE for object-oriented model-based reinforcement learning.
ACKNOWLEDGMENTS SA thanks to Kakao Brain and Center for Super Intelligence (CSI) for their support. ZL thanks to the ZJU-3DV group for its support. The authors would also like to thank Chang Chen and Zhuo Zhi for their insightful discussions and help in generating the 3D-Room dataset. A ADDITIONAL RESULTS OF SPACE

B ELBO DERIVATIONS
In this section, we derive the ELBO for the log-likelihood log p(x).
KL Divergence for the Foreground Latents Under the SPACE 's approximate inference, the D KL (q(z fg i |x) p(z fg i )) inside the expectation can be evaluated as follows.
KL Divergence for the Background Latents Under our GENESIS-like modeling of inference for the background latents, the KL term inside the expectation for the background is evaluated as follows.
Relaxed treatment of z pres In our implementation, we model the Bernoulli random variable z pres i using the Gumbel-Softmax distribution (Jang et al., 2016). We use the relaxed value of z pres in the entire training and use hard samples only for the visualizations.

C BOUNDARY LOSS
In this section we elaborate on the implementation details of the boundary loss. We construct a kernel of the size of the glimpse, gs × gs (we use gs = 32) with a boundary gap of b = 6 having negative uniform weights inside the boundary and a zero weight in the region between the boundary and the glimpse. This ensures that the model is penalized when the object is outside the boundary. This kernel is first mapped onto the global space via STN Jaderberg et al. (2015) to obtain the global kernel. This is then multiplied element-wise with global object mask α to obtain the boundary loss map. The objective of the loss is to minimize the mean of this boundary loss map. In addition to the ELBO, this loss is also back-propagated via RMSProp (Tieleman & Hinton. (2012)). This loss, due to the boundary constraint, enforces the bounding boxes to be less tight and results in lower average precision, so we disable the loss and optimize only the ELBO after the model has converged well. D IMPLEMENTATION DETAILS D.1 ALGORITHMS Algorithm 1 and Algorithm 3 present SPACE's inference for foreground and background. Algorithm 2 describe the details of rescale i function in Algorithm 1 that transforms local shift z shift i to global shiftẑ shift i .Algorithm 4 show the details of the generation process of the background module. For foreground generation, we simply sample the latent variables from the priors instead of conditioning on the input. Note that, for convenience the algorithms for the foreground module and background module are presented with for loops, but inference for all variables of the foreground module are implemented as parallel convolution operations and most operations of the background module (barring the LSTM module) are parallel as well.

D.2 TRAINING REGIME AND HYPERPARAMETERS
For all experiments we use an image size of 128 × 128 and a batch size of 12 to 16 depending on memory usage. For the foreground module, we use the RMSProp (Tieleman & Hinton. (2012)) optimizer with a learning rate of 1 × 10 −5 except for Figure 5, for which we use a learning rate of 1 × 10 −4 as SPAIR. For the background module, we use the Adam (Kingma & Ba (2014)) optimizer with a learning rate of 1 × 10 −3 . We use gradient clipping with a maximum norm of 1.0. For quantitative results, SPACE is trained up to 160000 steps. For Atari games, we find it beneficial to set α to be fixed for the first several thousand steps, and vary the actual value and number of steps for different games. This allows both the foreground as well as the background module to learn in the early stage of training.
We list out our hyperparameters for 3D large dataset and joint training for 10 static Atari games below. Hyperparameters for other experiments are similar, but are finetuned for each dataset individually. In the tables below, (m → n) : (p → q) denotes annealing the hyperparameter value from m to n, starting from step p until step q.

D.4 BASELINES
Here we give out the details of the background decoder in training of SPAIR (both full image as well as patch-wise training). The foreground and background image encoder is same as that of SPACE with the only difference that the inferred latents are conditioned on previous cells' latents as described in Section 2.2. For the background image encoder, we add an additional linear layer so that the encoded background latent is one dimensional.

SPAIR Background Decoder
Layer Size/Ch. Stride Norm./Act. For IODINE, we use our own implementation following the details as described in (Greff et al., 2019). For GENESIS, we also use our own implementation following the same architecture as in (Engelcke et al., 2019), but the details of individual networks are similar to SPACE's background module.

E DATASET DETAILS
Atari. For each game, we sample 60,000 random images from a pretrained agent (Wu et al., 2016). We split the images into 50,000 for the training set, 5,000 for the validation set, and 5,000 for the testing set. Each image is preprocessed into a size of 128 × 128 pixels with BGR color channels. We present the results for the following games: Space Invaders, Air Raid, River Raid, Montezuma's Revenge.
We also train our model on a dataset of 10 games jointly, where we have 8,000 training images, 1,000 validation images, and 1,000 testing images for each game. We use the following games: Asterix, Atlantis, Carnival, Double Dunk, Kangaroo, Montezuma Revenge, Pacman, Pooyan, Qbert, Space Invaders.
Room 3D. We use MuJoCo (Todorov et al., 2012) to generate this dataset. Each image consists of a walled enclosure with a random number of objects on the floor. The possible objects are randomly sized spheres, cubes, and cylinders. The small 3D-Room dataset has 4-8 objects and the large 3D-Room dataset has 18-24 objects. The color of the objects are randomly chosen from 8 different colors and the colors of the background (wall, ground, sky) are chosen randomly from 5 different colors. The angle of the camera is also selected randomly. We use a training set of 63,000 images, a validation set of 7,000 images, and a test set of 7,000 images. We use a 2-D projection from the camera to determine the ground truth bounding boxes of the objects so that we can report the average precision of the different models.


extfieldsofstudy:
Computer Science
Medicine

text:
Revealing the predictability of intrinsic structure in complex networks

Structure prediction is an important and widely studied problem in network science and machine learning, finding its applications in various fields. Despite the significant progress in prediction algorithms, the fundamental predictability of structures remains unclear, as networks’ complex underlying formation dynamics are usually unobserved or difficult to describe. As such, there has been a lack of theoretical guidance on the practical development of algorithms for their absolute performances. Here, for the first time, we find that the normalized shortest compression length of a network structure can directly assess the structure predictability. Specifically, shorter binary string length from compression leads to higher structure predictability. We also analytically derive the origin of this linear relationship in artificial random networks. In addition, our finding leads to analytical results quantifying maximum prediction accuracy, and allows the estimation of the network dataset potential values through the size of the compressed network data file.

P redicting structure or links of network is commonly defined as estimating the likelihood of existence of unobserved links or potential new links [1][2][3][4][5] . Network as a common form of data representation is ubiquitous across a broad range of fields from biology [6][7][8] , recommendation systems 9,10 to social media 1,3,11 . It represents the complex relationships or interactions among the elements in a system, which usually cannot be well described by simple mathematics. Hence, machine learning has been widely used in link predictions [10][11][12][13][14][15][16][17][18] . For instance, predicting protein-protein 7,8 or drug-target [19][20][21] interactions can guide more accurate biological experiments and reduce the experimental costs and time 2,6 . Despite the immense ongoing efforts in developing prediction algorithms, the fundamental understanding of the intrinsic prediction limit that can provide the much needed guidance is lacking. The difficulty lies with the fact that it is almost impossible to know the exact underlying mechanism of the network formation. In addition, the real networks are usually highly complex and gigantic in size, with many short feedback loops which cannot be effectively analyzed-a challenge 22 faced by both statistical physicists and computer scientists. Hence, understanding the prediction limits and quantifying it remains to be a long-standing challenge 2,23 .
In this study, we reveal the intrinsic predictability of networks through their structural properties. Intuitively, a network structure that can be captured with a few words means it is simple, and its links are easily predictable such as one-dimensional chain and two-dimensional lattice. Conversely, if a network requires lengthy description, then it has very complicated structure and its links are hard to predict. In computer language, the structure of any network can be encoded into binary strings. This motivates us to find the underlying relationship between the network compression length of the shortest binary string and the prediction limit. On the other hand, the inherent prediction limit is the maximum predictability or performance that the theoretical best predicting algorithm (TBPA) can achieve. Here the network structure predictability is defined as the performance of TBPA and is quantitatively measured using entropy in this paper. Thus, without knowing the exact underlying dynamics of network formation which determines this limit, we can use the best prediction algorithm available (BPAA) 1,2,23,24 to approximate the performance of TBPA. With these two quantities of the shortest lossless compression length and performance of TBPA, here, for the first time, we discover a linear relationship between them in different empirical networks such as biological networks, social networks, and technology networks (see Supplementary Note 1 for the detailed description of studied networks). Our finding implies that the shortest compression length of a network can tell us the structure predictability, which sets the limit of any prediction algorithm.

Results
Network shortest compression length. The shortest possible compression length can be calculated by a lossless compression algorithm 25 , which is a proven optimal compression for random networks and efficient for many real networks. Since all of the structural properties of two isomorphic networks are exactly the same, the algorithm first encodes the network structure into a string of binary codes, from which an isomorphic network can be reconstructed. To further remove the correlations in the structure, the binary string is compressed by a recursive arithmetic encoder 26 , which exploits the dependencies between the symbols in this binary string. After these two compression operations (see Fig. 1a), the length of the final bit string should be close to the network's structural entropy and can well measure its randomness (Shannon's source coding theorem 27 , see Supplementary Note 2). The length of this bit string is expected to increase as the structure becomes more random, and has been validated in Fig. 1b, which shows shuffled networks has longer compression length compared with the original empirical networks. As randomness in shuffling increases, the compression length increases monotonically.
Naturally the compression length is longer for a network with more nodes and links given the same randomness, i.e., the size of the network contributes to the compression length, rather than just the level of randomness of the network itself. In order to remove the size effect of the network, we normalize the compression length L through dividing it by the theoretical maximum compression length R, which corresponds to Erdős-Rényi (ER) 25,28 network of the same number of nodes and links (see Supplementary Note 2). The normalized value of the shortest compression length L * is given by Here R ¼ N 2 À Á hðqÞ À NlogN 25 , where N is the number of nodes and q is the probability having a link between any pair of nodes, also expressed as the binary entropy given by hðqÞ ¼ Àqlogq À ð1 À qÞlogð1 À qÞ.
Here log denotes the logarithmic operation with base 2, and we use this convention throughout this work.
Structure predictability measured by algorithm performance. Now we compare the shortest compression length with the performance of BPAA which approximates the structure predictability of network. In a large literature on link prediction algorithm [1][2][3][4][5][6]11,15,16 , it is common to assign each unlinked pair of nodes (possible missing link) a score, and higher score means higher likelihood of this pair being a missing link. Here we adopt the leave-one-out 31,32 approach and use the score of the existing links to quantify the BPAA performance. First we remove a link e i , and then use one particular prediction algorithm to estimate a score for each of the unlinked node pairs including e i . Based on the ranking of the scores in descending order, we get the ranking r i of the removed link e i . So when r i = 1, it means that the algorithm tells us the removed link is the most probable missing link among all other unlinked node pairs. We carry out this calculation of r i from the original network for every link at a time, and obtain a sequence of rank positions D = {r 1 , r 2 , . . . , r E }, where E is the total number of links in the original network. Naturally, the entropy of the distribution of D is a good holistic measure of the algorithm performance. For example, an ideal algorithm for a highly predictable network would have D = {r 1 = r 2 = ⋯ = r E = 1}, thus yielding the lowest distribution entropy of r i . Conversely, an ideal algorithm for a network with low predictability has very different values in D, leading to a high entropy of its distribution. In our calculation of this algorithm performance entropy H, the value r i can vary in the range 1 r i NðNÀ1Þ 2 À hkiN 2 þ 1 % N 2 2 , where NðNÀ1Þ 2 À hkiN 2 þ 1 is the total number of unlinked pairs, and 〈k〉 is the network's average degree. Thus we divide such range into bins with equal width N to avoid the contribution of network size N on the result (see Supplementary Note 3 for the discussion of other bin widths), and calculate H based on the probability distributions of the N∕2 bins: H ¼ À P N=2 j¼1 p j logp j where p j is the probability of r i in bin j (Supplementary Note 3). Figure 1c illustrates an example of such distribution for the Metabolic network 29 by RA algorithm 30 .
It is worth mentioning that using rank distribution entropy to measure the algorithm performance is the innovation in our work. In particular, we implement the leave-one-out method 31,32 to obtain the rank distribution. The other ranking methods 2,3,23 which pre-remove x fraction of links (x is usually 5% or 10%) with multiple samples are similar to the leave-one-out method, and would yield similar results as leave-one-out if x is very close to 0 (see Supplementary Note 4 and Supplementary Fig. 6). Using the leave-one-out method, there are three main advantages. Firstly, it is parameter-free, i.e., we do not need to consider the impact of the choice of x, and allows analytical study. Secondly, doing so preserves the original network structure as much as possible, such that the intrinsic true predictability is also preserved. Thirdly, the leave-one-out method has negligibly small fluctuation (see Supplementary Fig. 8). The main drawback is the higher computational complexity compared with only using 1 sample of removing x fraction of links. However, if multiple samples of removing x fractions are considered, the advantage of complexity of above other methods will not be significant over ours.
In order to obtain a good approximation of the TBPA performance of the network that is less dependent on a particular prediction algorithm being used, we employ a range of widely applied 11 prediction algorithms (see Supplementary Note 5), such as structural perturbation method (SPM) 23 , local random walk (LRW) 33 , average commute time 24 and common neighbor 1 to calculate the ranking positions r i , and use the one that gives the lowest H value H BPAA (Table 1) as the closest estimate of the network's predictability H TBPA . To further explore the relationship between network structure and H BPAA , it can be seen in Fig. 1d that as the network structure is progressively randomized due to shuffling, the entropy H BPAA increases, signifying a decrease in predictability. Eventually with enough shuffling, the algorithm performance entropy H approaches that of ER networks with a nearly uniform distribution of D (Fig. 1d inset). To further remove the effect of network size N and average degree 〈k〉, we normalize the BPAA performance entropy H BPAA value by logN À 1, which is the BPAA entropy of an ER network of the same size and average degree. Hence the normalized BPAA performance entropy H BPAA is defined as Empirical linear relationship. Having both network's structural shortest compression length and the BPAA performance, we are able to find the relationship between the two. Surprisingly, we discover a clear linear relationship between the H Ã BPAA and L * across 12 empirical networks in very different fields including biology, social media, transport and economics as shown in Fig. 2a: It is worth mentioning that, the way H and L are normalized is critical to obtain the linear relationship, as without proper normalization, the network's compression length cannot serve as a good indicator to quantify the structure predictability of a network (see Supplementary Note 6).
We find that such a linear relationship in Eq. (3) still exists even after we shuffle links as seen in Fig. 2b, or randomly add/ delete links as seen in Fig. 2c property of the shortest compression length 25,27 , i.e., the sum of the shortest compression lengths L A , L B of two independent networks A and B (with the same nodes but different links) is the same as the compression length L A+B of the combined network A + B (Fig. 2e). This is also valid for the compression algorithm in case of shuffling and adding links at random on empirical networks ( Fig. 2f, g). We now illustrate how to use the linear relationship in Eq. (3) to quantify the performance of a prediction algorithm. By comparing the Jaccard algorithm 35 performance entropy with the linear relationship (see Supplementary Fig. 11), it can be seen that for a given network, the further is the performance entropy from the linear line, the less optimal it is from the best possible prediction results (see Supplementary Note 7). This quantifiable distance can serve as a benchmark to practical algorithms in a given network. In other words, the further it is from the linear line, the more potential improvement it is possible for any new algorithm to achieve, warranting further effort on algorithm development. Conversely, if the algorithm performance entropy lies close to the line, it means that the algorithm is already performing quite well, and further exploration on algorithms is unlikely to yield significant improvement, suggesting diminishing returns in developing new algorithms.
Since the compression process can be executed through different node sequences, we also study the effect of compressed node sequences on the compression length, such as random, high degree priority and low degree priority. We find that the differences in compression lengths are very small (see Supplementary Fig. 2). Hence, we can see this compression algorithm is reliable when operating on real networks. Naturally, such relationship between H Ã BPAA and L * indicates that the link prediction limit of a network can be directly inferred from the network topological structures. More explicitly, by compressing the network structure into binary string and calculate the string's length, one can estimate the link prediction limit H Ã BPAA . Apart from the theoretical significance of the relation, one practical advantage of estimating H Ã BPAA with L * is that L * is only dependent on the network topological structure and has low computational complexity O(N + E) 25 . In contrast, estimation of H Ã BPAA requires a large number of available prediction algorithms, which usually have high computational complexity. Specifically, for calculating the rank distribution entropy H, we need to employ an algorithm for every link in a network to obtain its ranking. Consequently, the computational complexities of algorithms like SPM 23 and LRW 33 are about O(N 3 E) and O(N〈k〉 n E) in estimating H, where n is an integer constant. Furthermore, an additional benefit of the compression length is that it can be served as an independent indicator to identify missing or false links (see Supplementary Note 8).
Theoretical linear relationship. The empirical linear relationship inspires us to further figure out the underlying mathematical connection between the network shortest compression length and link prediction limit. For simplicity, we assume that an artificial network is generated from a static random matrix Q whose entry q ij denotes the link formation probability between node i and j. According to Shannon's source coding theorem 27 , the shortest compression length L of this artificial network is L ¼ P i>j hðq ij Þ À NlogN, which is the structural entropy 25 , where h(q ij ) is the binary entropy given by hðq ij Þ ¼ Àq ij logq ij À ð1 À q ij Þlogð1 À q ij Þ. To simplify representation, we use U to denote the term À P i > j q ij logq ij and expand the logarithm as a Taylor series, of which the higher-order small terms can be neglected, yielding: where ln denotes the natural logarithm, and we use this convention throughout this work. Since the occurrence of each link in this artificial network solely depends on the q ij , it is certainly true that the most accurate score that a TBPA could achieve should be exactly equal or proportional to q ij . Thus, without employing any prediction algorithm, we can obtain the ranking sequence D (see the right part of Fig. 3a) directly from Q and are able to numerically quantify H Ã TBPA in the same way that we calculate H Ã BPAA . We find that H Ã TBPA can be theoretically estimated by Q, given that the TBPA ranking distribution clearly agrees well with the distribution of q ij in Q (Fig. 3b). The entropy of the latter is given by Recalling that in the calculation of H BPAA , we have divided the ranking range into bins with equal width N. In such case, this coarse-grained distribution of Q is done by replacing every N values of q ij (in the descending order) by their average valueq ij (see the left part of Fig. 3a), and the  resulting entropy satisfiesH Q ¼ H Q À logN (see Supplementary Note 9). The difference between H TBPA andH Q is due to the fact that for H TBPA , link prediction algorithms only measure the likelihood of unobserved pairs of nodes forming a link. Such difference, however, contributes negligibly to the calculation of entropy as shown in Fig. 3b, yielding H TBPA %H Q . Taken together, we obtain (see details in Supplementary Note 9): We then normalize H TBPA by logN À 1 through the same way of Eq. (2), yielding Combining with Eqns. (1), (4)-(6) and eliminating the variable U, we obtain the linear relationship between L * and H Ã TBPA : When N is large Eq. (7) simplifies to In the thermodynamics limit Eq. (8) can be further approximated when logN ) loghki: Usually for many real networks, loghki logN is not negligible, and Eq. (8) is a better approximation. The detailed mathematics above is given in Supplementary Note 9.
To validate Eq. (7), we first construct artificial networks based on each empirical network in Fig. 2 with the same degree sequence, but a simple network formation mechanism based that Matrix Q distribution 5 6 TPBA performance entropy, H TBPA 7 8 9 Fig. 3 Theoretical relationship between the shortest compression length and network prediction limit. Networks are generated from a static random matrix Q with elements generated from the degree distribution of empirical networks. a Comparison between the ranking distribution from TBPA and that from random matrix Q based on Metabolic network 29 . The left half represents the ranking distribution of coarse-grained probabilitiesq ij , each of which is the average value of every N values of q ij ranked in a descending order. The right half illustrates the TBPA ranking distribution of existing links for this artificially generated network. The probability p j is an average value from 100 simulations. b Values of coarse-grained entropyH Q from the probability distributionq ij vs. TBPA performance entropy H TBPA for all of the artificial networks generated from the empirical ones. Each dot corresponds to the value of (H TBPA ;H Q ) of an artificial network, showing a good match between the two values. c Theoretical linear relationship between L * and H Ã TBPA calculated based on the empirical networks' degree distributions. Points of the same shape and color correspond to an empirical network (bottom-left) and its different shuffled versions, i.e., having its links shuffled different number of times up to all of original links. Each point represents the value pair of ðL Ã ; H Ã TBPA Þ of an artificial network with the same degree distribution as a corresponding (shuffled) empirical network. L * is calculated from Eq. (1) and H Ã TBPA is calculated from Eq. (6). The blue solid line is the average of 12 studied networks' analytical slopes obtained from Eq. (7) and the shaded region denotes the standard deviation. The gray dash line is the analytical result of Eq. (9) for the limit of 〈k〉 → ∞. d Competition between two terms of the slope in Eq. (8). e The empirical and theoretical value of slope in the linear relationship between L * and H Ã TBPA . The purple plane denotes the plane of empirical value at 1.63, and each colored point on the plane represents a real network. The lower curved surface below represents the theoretical values of the slope given by Eq. (7), and each colored point on the lower surface represents an artificial network constructed with the same degree distribution as original empirical networks, and the subsequent points with the same color are the artificial networks with randomly added links to increase 〈k〉. The red curve (lnN ¼ hki) on the lower surface is an estimated boundary (see Supplementary Note 9) that our theories are not valid far from the left side of it.
only depends on the static probability distribution matrix Q, with each element q ij ¼ k i k j 2E , where k i denotes the degree of node i in an empirical network. Secondly, we shuffle each above network by randomly picking a fraction of links and rewiring them randomly (same operation in Fig. 2b) and then construct artificial networks based on each shuffled empirical network's degree sequence. In the shuffling process, the randomness of the network increases with the increasing of the shuffling strength. Theoretically, the values (H Ã TBPA , L * ) changes along the straight line predicted by Eq. (7). It can be seen from Fig. 3c that this theoretical relation in Eq. (7) is well captured by our simulation for the different networks and their shuffled versions. We have also considered two other genuine edge-independent synthetic networks: the degree-correlated stochastic block model 36 characterizing the community structure, and the latent-geometric network model 37 characterizing network spatial structure. Both models have been shown to reproduce certain structural properties of real network. We find that the structure predictability properties obtained by these two models also follow our analytical result for artificial networks (see Supplementary Note 9 and Supplementary Fig. 15).
Surprisingly, the artificial networks generated by matrix Q leads to ðL Ã ; H Ã BPAA Þ pairs falling on a straight line as seen in Fig. 3c, with a slope that is different from the one given by Eq. (9) from the thermodynamics limit approximation. A careful examination shows that the artificial networks' average degrees 〈k〉 increase with the network sizes N faster than logN, balancing out the changes in both loghki logN and 2 hki in Eq. (8) (see Fig. 3d), therefore keeping the slope relatively constant around 0.7 instead of the thermodynamics limit 1 (see the plateau in Fig. 3e). For the approximation result in Eq. (9), we find that the thermodynamic limit slope 1 with large 〈k〉 can be realized in the artificial network with size larger than 10 100 (see Supplementary Fig. 14). But this large number is almost impossible in real networks. Moreover, the values of the slope between H Ã TBPA and L * of artificial networks are significantly lower than the empirical slope of 1.63 (Fig. 3e). This may be due to the more complex mechanisms and constraints in empirical network formations compared with simplistic random link connection in the artificial networks. Intuitively one expects the structure predictability of the empirical network to be higher than that of the purely random due to such additional constraints or mechanisms, and it is reflected in higher slope observed. Such hypothesis does not directly address the differences quantitatively, yet it sheds some light on the complex relationship between network entropy and structure predictability. At the same time, the empirical value of 1.63 hints at some strong universal mechanism that is common to all of these empirical networks.
Bounds of link prediction precision. In many practical applications involving link prediction, like recommendation system, a prediction algorithm gives the most likely missing links in the network among all possible links. One way to assess the algorithm's prediction precision is through the success of prediction [1][2][3][4][5][6]11,15,16 . In other words, we can look at the probability p 1 that the distribution of ranks given by the algorithm on the removed links defined in D. Note that if we take the top N predicted links to calculate the standard precision, then it is equivalent to p 1 used in this paper, since p 1 refers to the fraction of correct links among the top N predicted links. In practice, a good algorithm will give a ranking distribution that is usually monotonically decreasing, i.e., p 1 ≥ p 2 ≥ ⋯ ≥ p N∕2 ≥ 0 similar to Fig. 1b. That means the link predicted by the algorithm is indeed the most likely to be missing. Actually, we observe this phenomenon for all the algorithms and data used in this paper (see Supplementary Fig. 26). Using the linear relationship between L * and H Ã BPAA , for any network, we can calculate its normalized shortest compression length L * from its compressed binary string. Together with the assumption p 1 ≥ p 2 ≥ ⋯ ≥ p N∕2 ≥ 0, we arrive at an implicit function of the precision upper bound p 1 : The exact value of p 1 can be obtained by solving Eq. (10) (see Fig. 4a). Additionally, the lower bound p 1 of the prediction precision p 1 can be simply calculated through the following explicit formula (details in Methods section): One can define the prediction precision in a more general way other than p 1 . For instance, rather than defining the probability p 1 of the removed link falling into the first interval, one can loosen the definition to the first C intervals (1 ≤ C ≤ N∕2), i.e., defining precision as P C ¼ P C j¼1 p j (see Supplementary Note 10 for details on the derivations). Indeed, we validate the above upper bounds in the empirical networks as shown in Fig. 4b Commercial value of network dataset. In practice, the prediction bounds enable something that was not possible before, for instance it can be used to estimate the commercial value of a network dataset, through its compressed size without developing any prediction algorithm. Thus, for a compressed network with L bit, the conservative commercial value V is where Θ is an external economic variable. For example, in the scenario of inferring interactions between proteins, Θ can be written as θn, where θ represents the unit experimental costs that can be saved if successfully predicting one interaction and n(≤N) denotes the number of predicted interactions. In other words, once we obtain the smallest number of bits L needed to store a compressed network data file, we can directly derive the approximate value of this dataset through Eq. (12). Note that the commercial value we refer to above is the potential additional value that can be realized by predicting unobserved information from the network, without considering the overlapping value from external information outside the network structure. A more general framework of data commercial values when various external information is available is discussed in Supplementary Note 11.

Discussion
Although we have seen that our finding applies to a broad range of empirical networks and their artificial counterparts, it is important to take note of certain limitations in practice. It is know from ref. 38 that Eq. (4) is accurate when lnN ( hki ( N À lnN (see the red curve lnN ¼ hki in Fig. 3d). Therefore when a network is extremely sparse or dense, the relationship between L * and H Ã BPAA may not be close to what we found. In addition, the network entropy estimation and structure predictability analysis assume randomness in structure. That means for regular networks like lattice structure, our finding does not hold. To illustrate the impact of regular structural features, we combine a regular network from a circular model network 39 (Fig. 5a) and real networks. The analysis on such synthetic networks show that indeed, such regular structure in a random network makes the H Ã BPAA vs. L * relation deviates from the slope of 1.63 as seen in Fig. 5b. Hence our result is more valid for networks without significant amount of regular links. That being said, if a network is very regular in structure, lattice for instance, its structure is easy to predicted and compressed due to high regularity involved.
In conclusion, for the first time we established a theoretical framework to quantify the intrinsic structure predictability in real and artificial random networks solely based on network structure, and independent from any prediction algorithms. With theoretical intuition on proper normalization, we uncover a universal linear relationship between the shortest compression length of the network structure and its structure predictability for a wide range of real networks such as biological networks, social networks and infrastructure networks, and analytically derived the mathematical origin of this linear relationship in artificial networks. In principle, such relationship can serve as a benchmark to quantify the performance of any practical prediction algorithm in nonregular networks. Leveraging upon this linear relationship, we can obtain the structure predictability that is intrinsic to the structure of complex networks and provide the accuracy bounds for any link prediction algorithm. In practice, our method can also be used to estimate the commercial value of a network dataset through its compressed length without using any link prediction algorithm. Our finding is demonstrated upon data structure in the form of networks. However, if the structural entropy and prediction limit are linked through the underlying dynamical process for data structures other than networks, it is likely that the predictability in some other machine learning problems for different types of systems can also be inferred through similar approaches of optimal compression.

Methods
Compression algorithm. Our compression scheme builds upon the seminal paper of 25 , which is a two-step lossless compression of graphs. First we encode a network into two binary sequences B 1 , B 2 , through traversing all nodes on the network from an initial node, and encode nodes' neighbors according to specific rules 25 . In the second stage, both B 1 , B 2 are compressed by an improved arithmetic encoder based on recursive splitting proposed by K. Skretting 26 to exploit the dependencies between the symbols of the sequence. After that we obtain two compressed binary sequencesB 1 ,B 2 , and we define the compression length of the network L to be the total length of these two sequences: where 'ðB 1 Þ,'ðB 2 Þ are the length ofB 1 andB 2 , respectively. More details about the algorithm is provided in Supplementary Note 2.
Upper and Lower Bounds of p 1 and P C . Here we demonstrate how to identify the upper and lower bounds of link prediction precision p 1 and P C . The basic idea is to transform these problems into optimization problems of certain boundary conditions. Firstly, we calculate the network's BPAA performance entropy H Ã BPAA through its shortest compression length L * given by Eq. (3). Its un-normalized value is ð1:63L Ã À 0:63ÞðlogN À 1Þ, and by the definition of entropy it can be written as There are two constraints for p i s. The first is that they sum up to unity, i.e.   and: as they are arranged in decreasing order. With these boundary conditions, p 1 is maximized when the other probabilities are equal. Then it can be found that the upper bound of p 1 is given by Eq. (10).
And the minimum value of p 1 is when there are as many p i s as close to p 1 as possible, leading to a lower bound given by Eq. (11).
For predictability of the top C intervals P C ¼ P C j¼1 p j , the upper bound corresponds to the case when the last N∕2 − C probabilities are the same, yielding ÀP C log P C C À ð1 À P C Þ log 1 À P C N=2 À C ¼ ð1:63L Ã À 0:63ÞðlogN À 1Þ: For the lower bound of P C , the value of p C can be introduced to construct the minimization problem. It has the boundary condition of 0 ≤ p C ≤ P C ∕ C. The approximate lower bound is the solution of the following: P C % minfP C : ÀγðP C ; p C ÞlogγðP C ; p C Þ À ð1 À γðP C ; p C ÞÞlog p C ¼ ð1:63L Ã À 0:63ÞðlogN À 1Þg: The more detailed mathematics is provided in Supplementary Note 10.

Data availability
Data of this study is available at http://www.huyanqing.com/.