A Tensor-Based Framework for Studying Eigenvector Multicentrality in Multilayer Networks
Mincheng Wu Shibo He Yongtao Zhang Jiming Chen Youxian Sun Yang-Yu Liu Junshan Zhang H. Vincent Poor


Low Cost, Open-Source Testbed to Enable Full-Sized Automated Vehicle Research
    Austin Costley Chase Kunz Ryan Gerdes Rajniaknt Sharma 
    
=============================================================================




 
  * State Key Laboratory of Industrial Control Technology, Zhejiang University, Hangzhou 310027, Zhejiang, China.
 
  * Channing Division of Network Medicine, Brigham and Women's Hospital and Harvard Medical School, Boston, MA 02115.
 
  * Center for Cancer Systems Biology, Dana-Farber Cancer Institute, Boston, MA 02115.
 
  * School of Electrical, Computer and Energy Engineering, Arizona State University, Tempe, AZ 85287.
 
  * Department of Electrical Engineering, Princeton University, Princeton, NJ 08544.



()





Abstract



Centrality is widely recognized as one of the most critical measures to provide insight in the structure and function of complex networks.
While various centrality measures have been proposed for single-layer networks, a general framework for studying centrality in multilayer networks (i.e., multicentrality) is still lacking.
In this study, a tensor-based framework is introduced to study eigenvector multicentrality, which enables the quantification of the impact of interlayer influence on multicentrality, providing a systematic way to describe how multicentrality propagates across different layers.
This framework can leverage prior knowledge about the interplay among layers to better characterize multicentrality for varying scenarios.
Two interesting cases are presented to illustrate how to model multilayer influence by choosing appropriate functions of interlayer influence and design algorithms to calculate eigenvector multicentrality.
This framework is applied to analyze several empirical multilayer networks, and the results corroborate that it can quantify the influence among layers and multicentrality of nodes effectively.







§ INTRODUCTION

Centrality quantifies  the importance of nodes in a graph, and has been widely studied to understand the structure and function of complex networks<cit.>.
For example, it can be utilized to  identify the most influential person in an online social network<cit.>, the most crucial artery in transport congestion<cit.>, or the most important financial institution in the global economy<cit.>.
Over 30 different  centrality measures (e.g., degree centrality, betweenness centrality, closeness centrality,  eigenvector centrality, and control centrality) have been examined  in the literature<cit.>.
Among these, eigenvector centrality, defined as the leading eigenvector of the adjacency matrix of a graph, has received  increasing attention<cit.>.
It is worth noting that PageRank, a variant of  eigenvector centrality,
is the primary algorithm used in Google's search engine<cit.>.

Notably, most previous studies have focused on eigenvector centrality in a single-layer network, in which all nodes/links are assumed to be of the same type (centrality-homogeneous).
As revealed recently<cit.>, many practical complex systems, ranging from the Internet to airline networks, have multiple types  of nodes and/or  links between nodes.
Multilayer networks, which consist of multiple layers of nodes with intra- and interlayer links, can be used to model such  complex systems.
Figure <ref> shows two examples of multilayer networks (see Supplementary Fig. 1 for more examples).
Simply aggregating a multilayer network into a single-layer one would obviously lead to a miscalculation of centrality.
Recent work on eigenvector-like centrality in multilayer networks either  assigned constant weights to predetermine interlayer influence (which can be regarded as the gain or loss of the interplay strength between two layers<cit.>), or focused on a special case of multilayer networks, i.e., the so-called multiplex networks (where all layers share the same set of nodes, and  interlayer links only exist between counterpart nodes)<cit.>.
It is of significant interest to develop a framework for studying eigenvector-like centrality in general multilayer networks, hereafter  referred to as eigenvector multicentrality.
In this study, we introduce a tensor-based framework that enables the quantification of the relationship between interlayer influence and eigenvector multicentrality.
It is challenging to compute eigenvector multicentrality of nodes in such a framework since  interlayer influence  and eigenvector multicentrality are interdependent.
We prove the existence and uniqueness of eigenvector multicentrality for given appropriate forms of interlayer influence.
We also design efficient algorithms to calculate it for two interesting scenarios.
This framework offers a novel approach for  modeling and quantifying the interlayer interactions in multilayer networks, providing a systematic way of characterizing eigenvector multicentrality.
Experimental results based on several real-world multilayer networks corroborate our analytical results.



§ RESULTS



 §.§ A tensor-based framework for studying eigenvector multicentrality.


In the calculation of eigenvector centrality of nodes in a single-layer network, a directed link to a node can be viewed as  a vote of support.
Each node fairly propagates its entire centrality score to its neighbors recursively. 
The eigenvector centrality of a node is defined as the scores that it gathers  from its neighbors after appropriate normalization in the steady state.
Formally, the vector consisting of the eigenvector centrality of all nodes is defined as the leading left eigenvector of the adjacency matrix associated with the single-layer network.
We generalize this definition to multilayer networks by taking into account interlayer influence among layers.
Specifically, a multilayer network is modeled as ℳ=(ℒ,ℰ), where ℒ={L_α;α=1,2,⋯, K} is a collection of graphs L_α=(V_α, E_α) representing layers in ℳ;  V_α={v_1,α,v_2,α,⋯ v_n_α,α} is the set of nodes, where n_α denotes the number of nodes, and E_α is the set of intralayer links in layer α; ℰ={E_αβ⊆ V_α× V_β; α,β =1,2,⋯,K (αβ) } contains the interlayer links in ℳ<cit.>.
To avoid confusion, we use Latin letters {i,j,⋯} to indicate nodes and Greek letters {α,β,⋯} to indicate layers.
Tensors provide a general mathematical tool to describe high-dimensional objects<cit.>, and notably, tensors have been employed to study multilayer networks<cit.>.
For example, a fourth-order tensor M_jβ^iα, called the adjacency tensor, is used to encode a directed, weighted link from node i in layer α to node j in layer β (see Supplementary Text I for further details on tensorial representations).
We further introduce the influence tensor W^α_β, which is a second-order tensor measuring the interlayer influence from layer α to layer β.
In our framework, the influence tensor W may be treated as a constant tensor when quantitative knowledge is available.
The interaction tensor is defined as H^iα_jβ=W^α_β M_jβ^iα, encoding the interaction from node i in layer α to node j in layer β.
Intuitively, when the influence from  layer α to layer β is greater than one, the centrality scores propagating along the links from layer α to  layer β will be magnified,  and vice versa.


The second-order tensor Φ_iα is defined as the solution to the following tensorial equation:

    H_jβ^iαΦ_iα=λ_βΦ_jβ,

where λ_β is a coefficient related to layer β,
and the Einstein notation<cit.> is adopted here (see Supplementary Text I for further details).
Because Φ is a eigenvector-like centrality, it will hereafter be referred to  as eigenvector multicentrality, and Φ_iα represents the eigenvector multicentrality score of node i in layer α.
In  the calculation of  eigenvector multicentrality, after each node propagates  its entire multicentrality score to neighbors, the scores from layer α to layer β will be multiplied by the influence coefficient W^α_β. 
Hence, we can obtain the eigenvector multicentrality via appropriate normalization in the steady state.
Notice that the normalizing coefficient λ_β may be different for different layers in a multilayer network. 
This differs from the eigenvector centrality in a single-layer network, where all nodes share a common normalizing coefficient λ_1 (namely, the leading eigenvalue of the adjacency matrix).



Many existing models can be incorporated in our eigenvector multicentrality framework  by choosing the influence tensor W appropriately.
For instance, in an author-document heterogeneous network, unweighted interlayer links connect documents to their authors.
A directed unweighted intralayer link exists between two documents if one document refers to the other and the undirected weighted intralayer link between two authors represents their social tie.
The multicentrality in such an author-document network was defined as the leading eigenvector of a stochastic matrix, which encodes the probabilities that a random surfer moves along intra- and interlayer links  in a combined random walk process<cit.>. Clearly, this model can be incorporated in our framework by setting the intra- and interlayer influence weights to predetermined constants.
Further, a definition of multicentrality in multiplex networks has been proposed in Ref. 20, considering the influence from counterpart nodes in other layers by importing the influence matrix Q=(q_αβ)∈ℝ^K× K. The influence matrix is non-negative, and q_αβ measures the influence of layer β on layer α.
One can calculate eigenvector-like multicentrality of a multiplex network once Q has been obtained. 
Observe that in Ref. 20 the influence matrix Q is predetermined and given as in a lower-dimensional form (i.e., a matrix), whereas in our framework the influence tensor W depends on the eigenvector multicentrality and hence they are interdependent.
Moreover, interconnected multilayer networks have been proposed to predict diffusive and congestion processes<cit.>, where the undirected unweighted interlayer links connect nodes to their counterparts in other layers. 
Eigenvector centrality in these networks  is a special case of our framework where the interlayer influence is equal to one.



 §.§ Leveraging prior knowledge about interlayer interactions.

Quantifying the influence tensor W is important to calculate multicentrality in our framework.
As expected, the influence tensor W is typically a function of the adjacency tensor M and the multicentrality Φ, rather than being a constant.
In practice, precisely predetermining the influence tensor is often infeasible.
One advantage of our general framework lies in the leveraging of prior knowledge about the influence tensor W in diverse applications and the  calculation of multicentrality even when W and Φ are interdependent.

Consider a typical scenario in which all nodes are centrality-homogeneous and layers are heterogeneous in a multilayer network.
In such a scenario, the multicentrality scores of all nodes are comparable, and can be represented by a vector C∈ℝ^N, where N is the number of all nodes. We call the eigenvector multicentrality in such a scenario  global multicentrality. 
We normalize the multicentrality such that the vector C is defined over an N-dimensional simplex.
For example, in a multilayer network consisting of web pages on different subjects (see the empirical results for further details), we may need to compare the multicentrality scores of two web pages on different subjects.
Clearly, the multicentrality score propagates differently along interlayer links and along intralayer links, owing to differences in the popularity of different subjects. 
Here, we assume that the importance of layer α is a function of the multicentrality of all nodes in layer α, denoted by f(Φ_:α), where the colon “:" indicates all elements of a given dimension<cit.>. 
Notably, the function f(·) describing the layer importance is application-dependent.
The function f could be, for example, the L^1-norm f(Φ_:α)=Φ_:α_1, which means the aggregated multicentrality scores of nodes in layer α, or it could be f(Φ_:α)=Φ_:α_1/n_α, which denotes the average multicentrality score over nodes in layer α.


There are a variety of ways to define the influence tensor. In this article, we define

    W^α_β=f(Φ_:α)/f(Φ_:β).

That is to say, the interlayer influence between two layers depends on their relative layer importance.
It is clear that there is  no gain or loss for links between nodes in the same layer,  because the interlayer influence W^α_α (α=1,2,⋯,K) is equal to one.
For  links from a node in a more important layer, there is a gain in the multicentrality, and vice versa.
Then the interaction tensor can be written as  H^iα_jβ=f(Φ_:α)/f(Φ_:β)· M_jβ^iα. Further, we prove that λ_α=λ_1, ∀α∈{1,2,⋯, K}, and Eq. <ref> reduces to

    H_jβ^iαΦ_iα=λ_1Φ_jβ,

where λ_1 is the leading eigenvalue of the interaction tensor H, and is irrelevant to β in this scenario.
Notably, λ_1 is also the leading eigenvalue of the adjacency tensor M, indicating that the interaction tensor H maintains the leading eigenvalue of the adjacency tensor M (see Supplementary Text V for further details).
From Eq. <ref>, we can see how multicentrality score propagates in a multilayer network.
Considering that a node i in layer α links to another node j in layer β, node i will propagate its multicentrality score to node j scaled by an influence coefficient W^α_β, which could be a gain (W^α_β>1), a loss (W^α_β<1), or even (W^α_β=1), and the multicentrality Φ_jβ comprises the scores that node j in layer β gathers in the steady state.
The existence and uniqueness of the multicentrality Φ are proved in Supplementary Text V. 
For a given function f, we can calculate the global multicentrality using the compressed power iteration method introduced in the Materials and Methods section.


We consider another interesting scenario in which nodes in different layers are  heterogeneous, and thus are not comparable. For example, in a heterogeneous network consisting of authors and papers, it is not meaningful to compare the multicentrality of an author with that of a paper.
Because the multicentrality of nodes in different layers may have varying implications, we can only calculate the local multicentrality of nodes in each layer while taking into account the interlayer influence, where local multicentrality means that the nodes in each (local) layer are centrality-homogeneous.
In such a scenario, the multicentrality score of a node cannot simply propagate along interlayer links to other layers.
We measure the local multicentrality of nodes in each layer by defining the influence tensor  in the framework as

    W^α_β=∑_i,j=1^NM^jβ_iαΦ_jβ/∑_i,j=1^NM^iα_jβΦ_iα.

Notice that the denominator ∑_i,j=1^NM^iα_jβΦ_iα in W^α_β is a normalizing constant, which is the sum of scores propagating along interlayer links from layer α to layer β. 
Moreover, the sum of interactions from nodes in layer α to nodes in layer β is given by the numerator ∑_i,j=1^NM^jβ_iαΦ_jβ, which is the sum of scores propagating from layer β to layer α.
Therefore, by defining the influence tensor in Eq. <ref>, we assume that the score flow going out of one layer will be returned to the layer.
In such a way, we can calculate the local multicentrality of nodes in each layer independently while the interlayer influence is taken into account.
The detailed proofs of the existence and uniqueness of Φ in local multicentrality are provided in Supplementary Text VI.
We can also calculate the local multicentrality Φ numerically using the compressed power iteration method.


Because the prior knowledge is network-specific and application-dependent, we present two interesting scenarios, for global multicentrality and local multicentrality respectively, to illustrate how to leverage prior knowledge to find W and compute the multicentrality of nodes.
The PageRank algorithm has been now widely utilized in social, transportation, biology and information network analysis for link prediction, recommendation, etc.<cit.>.
Note that PageRank centrality is a variant of eigenvector centrality.
In PageRank centrality, each node distributes its PageRank score to its neighbors along outgoing links on an equal footing, and a node's PageRank score is defined as the sum of scores that it gathers from its neighbors in the steady state.
Our eigenvector multicentrality framework can be easily carried over to characterize PageRank multicentrality (see Supplementary Text IV for further details about PageRank multicentrality).




 §.§ Multicentrality in empirical networks.


We first  consider a dataset from Wikipedia consisting of 4,604 web pages (see the Materials and Methods section for further details about this dataset and how it has been obtained).
The web pages are divided by Wikipedia into 15 subjects, including art, business studies, citizenship, countries, design and technology, everyday life, geography, history, IT, Language and literature, mathematics, music, people, religion and science, and one web page belongs solely to one subject.
We build a multilayer network by placing web pages of the same subject in the same layer, and establishing a directed link between two web pages if there is a hyperlink between them (a subnetwork is shown in Fig. 1A).
For comparison, we also build a single-layer network by aggregating all web pages in all layers (see Supplementary Fig. 4 for further details).
We select 4,128 web pages from these, in order to guarantee that the network is connected and all nodes have at least one out-degree and one in-degree.

We measure the PageRank multicentrality of web pages in the constructed multilayer network.
Moreover, because web pages in different subjects are centrality-homogeneous, global PageRank multicentrality is used.
We consider three common forms for layer importance f(·):
f_1(Φ_:α)=ln(1+N·|Φ_:α|_1/n_α), f_2(Φ_:α)=|Φ_:α|_∞ and f_3(Φ_:α)=|Φ_:α|_1/n_α.
Table 1 shows the results of global PageRank multicentrality (f=f_1) in the Wikipedia multilayer network and PageRank centrality in the aggregated Wikipedia network, where the digits in parenthesis indicate the differences between these two rankings (see Supplementary Tables 1-3 for further details).

Note that the entry “United States” has the largest multicentrality score, because  it has the largest number of incoming links. Furthermore, the  web pages linked to it have high multicentrality scores.
The entry “Europe” in the layer “Geography” has many enhanced links from  nodes in the layer “Countries”, because the average PageRank multicentrality score of nodes in “Countries” is higher than that in “Geography”.
“France” has  large numbers of  incoming links from entries in the layers “Art”, “Music” and other subjects; scores along these incoming links, however, will be diminished since “Art” and “Music” are less important than “Countries”.
Note that the entry “Television” is significantly promoted, because the layer importance of “Design and Technology” is relatively low. 
Global multicentrality can effectively quantify the influence between layers even when we have limited prior knowledge, rather than aggregating all nodes without considering interlayer influence, as shown in many previous methods.

Eigenvector multicentrality is an effective predictor for searching for the most important node in multilayer networks, which is also validated by the results obtained from the page views (PVs) in Wikispeedia.
Wikispeedia is a human-computation game<cit.>, in which users are requested to navigate from a given web page to a  target one by only clicking on Wikipedia links. 
We collect all completed navigation paths and obtain the PVs of each web page from Wikispeedia.
For all entries, we calculate their PageRank multicentrality scores, PageRank centrality and degree centrality scores, and compare them to their PVs  in Wikispeedia (see Supplementary Tables 4-8 for further details). 
The results are shown in Figs. <ref> and <ref>.
We also list the average PageRank multicentrality score of each layer (subjects) in Wikipedia (see Supplementary Tables 9-11 for further details).
The Spearman rank correlation  coefficients show that PageRank multicentrality outperforms PageRank centrality and degree centrality in the aggregated network.

Next, we consider a transportation network consisting of airports and air routes between them. We first consider 450 airports in Europe<cit.> (see the Materials and Methods section for details about this dataset and how it has been obtained), and we focus on three main airlines. 
We then build a multiplex network with three layers (airlines) and  450 nodes (airports) in each layer as shown in Fig. <ref>B, where the dotted lines are interlayer links between airports and their counterparts in other layers.
We measure the global PageRank multicentrality in the multiplex network using three forms of layer importance: f_1(Φ_:α)=e^|Φ_:α|_1/n_α-1, f_2(Φ_:α)=|Φ_:α|_1/n_α
and f_3(Φ_:α)=ln(1+N·|Φ_:α|_1/n_α).
Then, the PageRank multicentrality score of each airport is obtained by assembling the multicentrality scores of all its counterparts in all layers.
For comparison, we also build a single-layer network, called the aggregated network, by combining the same airports in the three layers.
We further introduce the versatility, a good predictor for diffusive and congestion processes in multilayer networks<cit.>, which is a special case of our framework when setting  all components  of the influence tensor W to one.
We focus on the coverage ρ(t), a suitable proxy for the exploration efficiency of the network<cit.>, defined as the average fraction of distinct nodes being visited up to time t regardless of the layer,  assuming that a walker starts from a certain node in the network
(see the Materials and Methods section for further details about the coverage).
We investigate whether  multicentrality helps understand the role that a node plays in dynamical scenarios.
To this end, we compute the Spearman correlation coefficient between the ranking of the airports by multicentrality and that by the coverage at time t of a hypothetical epidemic spreading process that starts from a certain airport.
For comparison, we also compute baselines such as the rankings by versatility and PageRank centrality in the aggregated network (see Supplementary Tables 12-14 for further details).
We calculate the Spearman correlation coefficients for these five methods at each time step, and the results are shown in Fig. <ref>, where the time ranges from t=1 to t=4,000.
It is shown that the three multicentrality measures achieve higher accuracy (their  correlation coefficients exceed 0.947) in the steady state (t≥3,000).
We then perform a similar analysis on an airline network from the U.S., which contains the airlines flying from the U.S. on Jan 3, 2008 (with data provided by the American Statistical Association Sections on Statistical Computing, <http://stat-computing.org/>). 
We build a multiplex network with 20 layers and  284 nodes in each layer, and a similar conclusion can be drawn (see Supplementary Table 15 for further details).

Another real-world example we consider is a social network, where we apply our framework to an  e-mail network constructed from a large European research institution with 1,005 nodes (individuals) and 42 layers (departments), on which we consider an epidemic spreading process. 
The simulation results indicate that the nodes with higher eigenvector multicentrality play a more important role in the  epidemic spreading process (see Supplementary Text VII for details).



§ DISCUSSION

As shown in recent work on eigenvector centrality (and its variants)<cit.>, it is of significant interest to build a framework for studying eigenvector-like centrality in multilayer networks.
The existing studies, however, assumed empirical influence coefficients or relied on specific types of multilayer networks.
Here we develop a general framework for studying eigenvector multicentrality in multilayer networks, which enables the quantification of the impact of interlayer influence on eigenvector multicentrality, providing an analytical tool to describe how eigenvector multicentrality propagates among different layers.
Further, this framework can easily leverage prior knowledge about the interplay among layers to characterize eigenvector multicentrality for varying scenarios.
As the interlayer influence and multicentrality of nodes are interdependent, they are jointly solved using a compressed power iteration method.
Furthermore, we formulate and analyze PageRank multicentrality for practical applications within the proposed framework.
We also perform theoretical analyses to prove the existence and uniqueness of the solutions in Supplementary Text V-VI, which allows us to calculate global and local multicentrality in any strongly connected multilayer networks. 
For an arbitrary multilayer network, we treat a dead end (a node with no outgoing links) the same as if it had outgoing links to all nodes, and introduce a damping factor to guarantee the existence and uniqueness of the solution.
The results from  empirical networks demonstrate that our general framework can effectively quantify the interlayer influence, and eigenvector multicentrality is a good measure to identify important nodes from both structural and dynamical perspectives.
Thus,  multicentrality aids in understanding and predicting the behaviors of dynamic processes by leveraging network structure, and describes the structure-function relationship of multilayer networks well.


We believe that the concept of multicentrality has the potential to enable a deep understanding of the structure and function of multilayer networks.
Because the real-world scenarios of multilayer networks vary, a key step is to find appropriate forms of interlayer influence for theoretical analyses. Here we  consider two interesting scenarios, for global multicentrality and local multicentrality respectively.
We believe that our proposed tensor-based framework can be applied to  more empirical networks in various  scenarios, including social networks, transportation networks, biological networks, etc.




§ MATERIALS AND METHODS



 §.§ Numerical solution.

The crux of the proposed framework is to solve the tensorial equation

    H_jβ^iαΦ_iα=λ_βΦ_jβ,

where the solution Φ is  the multicentrality tensor.
To obtain the numerical solution, we first flatten the adjacency tensor M into a matrix, i.e., we represent the fourth-order tensor M∈ℝ^N× N× K× K as a matrix M∈ℝ^NK× NK, where M denotes the lower-dimensional form of the tensor M.
Then we vectorize the second-order tensor Φ∈ℝ^N× K into a supravector Φ∈ℝ^NK and denote by W the matrix form of the influence tensor W (see Supplementary Text II for further details of  tensor decomposition). Further, we denote by Λ=[λ_1,λ_2,⋯,λ_K]^⊤∈ℝ^K  a vector encoding the  normalizing coefficient in each layer.
Thus, we obtain the matrix equation

    (W⊙M)^⊤·Φ=Λ⊙Φ,

where ⊙ denotes the Khatri-Rao product and W=W(M,Φ) is  a function of the matrix M and the multicentrality Φ.


For the numerical solution, we propose a compressed power iteration method, whose iteration scheme is as follows:

    Φ^(k+1)=Φ^(k)+D^(k)[(W⊙M)^⊤(k)⊙Ω^(k)-E_N]Φ^(k),

where Ω^(k)=[ω^(k)_1,ω^(k)_2,⋯,ω^(k)_K]∈ℝ^K is the normalizing vector. Denoting B^(k)=(W⊙M)^⊤(k)⊙Ω^(k)-E_N, we can write the iteration scheme as

    Φ^(k+1)=Φ^(k)+D^(k)B^(k)Φ^(k),

where D^(k)∈ℝ^N× N is a diagonal matrix related to B^(k), and D^(k) compresses the induced infinity norm of the matrix B^(k) such that the L^1-norm of each row in D^(k)B^(k) is strictly less than one.

Specifically, we let ω^(k)_γ=Φ_γ^(k)_1^-1. With regard to the global multicentrality, the vector Λ contains equivalent elements, i.e., Λ=[λ_1,λ_1,⋯,λ_1]^⊤. Thus we have Ω^(k)=[ω^(k)_1,ω^(k)_1,⋯,ω^(k)_1]^⊤, where ω^(k)_1=Φ^(k)_1^-1.
Further, we can specify the diagonal matrix D^(k) to compress the induced infinity norm of matrix B^(k), where  D^(k) is not unique in practice. For example, we could take  D^(k)=( diag{(B^(k)+E_N)·1_N× 1})^-1.
Then for each Φ^(k)>0 (i.e., all the components in Φ^(k) are positive), the matrix [E_N+D^(k)B^(k)] is strictly diagonal dominant with positive elements.
Hence, the iterations in Eq. <ref> converges to a unique solution<cit.> and this solution satisfies the matrix Eq. <ref> (see Supplementary Text VIII for further details).



 §.§ Multilayer network of  Wikipedia.

The Wikipedia dataset  contains 4,604 entries and 119,882 hyperlinks<cit.> (data provided by the Stanford Network Analysis Project, <http://snap.stanford.edu/index.html>).
To ensure connectivity, we select the entries that have at least one out-degree and one in-degree. Then 4,128 entries and 113,441 links are obtained and these entries are divided into 15 subjects by Wikipedia.
Resultantly, we can build a multilayer network with N=4,128 nodes and K=15 layers, where each layer contains the entries of a single subject. The adjacency tensor M ∈ℝ^N× N× K× K encodes the directed and unweighted links of the multilayer network and the influence tensor W ∈ℝ^K× K encodes the interlayer influence between any two layers.
We measure the interlayer influence

    W^α_β=f(Φ_:α)/f(Φ_:β),

where f(Φ_:α) indicates the layer importance of layer α.
Here, we consider three forms of layer importance:
f_1(Φ_:α)=ln(1+N·|Φ_:α|_1/n_α), f_2(Φ_:α)=|Φ_:α|_∞ and f_3(Φ_:α)=|Φ_:α|_1/n_α.
Further, we obtain the interaction tensor H ∈ℝ^N× N× K× K as H^iα_jβ=W^α_β M_jβ^iα.
Following the construction of the interaction tensor H, we then solve the tensorial equation

    H_jβ^iαΦ_iα=λ_1Φ_jβ

using the compressed iteration method.
Finally, the multicentrality tensor Φ is in the space ℝ^N× K, and  Φ_iα represents the multicentrality of  node i in  layer α.



 §.§ Multilayer network of the European airlines.

The European airline network contains  450 airports in Europe and the air routes for 37 airlines (see Ref. 34 for more details about this dataset).
For each airline, we can build a network with N=450 nodes and a set of links representing routes between  airports.
We select those airlines with the number of air routes greater than N/2, such that the average degree for each node  in the constructed network is at least one. 
In this way, we obtain three main airlines: Ryanair, Lufthansa and Easyjet.
We thus have a multiplex network with K=3 layers and 450 nodes in each layer.
Then, we interconnect  the same airport across layers, obtaining a three-layer multiplex network.
The adjacency tensor M is in the space ℝ^N× N× K× K.
In the context of the  multiplex network, M_jα^iα encodes the undirected and unweighted intralayer links in layer α, while M_jβ^jα=1 encodes the undirected and unweighted interlayer link for node j between layer α and layer β.
For the interlayer influence, we again consider three forms of the layer importance:
f_1(Φ_:α)=e^|Φ_:α|_1/n_α-1, f_2(Φ_:α)=|Φ_:α|_1/n_α
and f_3(Φ_:α)=ln(1+N·|Φ_:α|_1/n_α). After obtaining the influence tensor W via  Eq. <ref>, we have the interaction tensor H^iα_jβ=W^α_β M_jβ^iα in the space ℝ^N× N× K× K.
Finally, we solve the tensorial Eq. <ref> using the compressed iteration method and obtain the multicentrality tensor Φ_iα.

With respect to the coverage ρ(t), we have

    ρ(t)=1-1/N^2∑^N_i,j=1δ_i,j(0)exp[-P_j(0)ℙE_i^⊤],

where δ_i,j(0)=0 for j=i, and δ_i,j(0)=1 otherwise. Here, P_j(0) represents the supravector of probabilities at time t=0 (assuming that the walker starts at node j) and the matrix ℙ indicates the probability of reaching each node through any path of length 1,2,⋯,t+1.
Furthermore, E_i=(e_i,e_i,⋯,e_i) is the supravector in which e_i is the ith canonical row vector (see Ref. 35 for details about the derivation of Eq. <ref>).

























§ SUPPLEMENTARY INFORMATION



 §.§ Supplementary Text



  §.§.§ I. Tensorial representations

We first introduce a mathematical formulation of multilayer networks<cit.>, in which tensorial representations are adopted.
The column vector e(s)=(0,…,0,1,0,…,0)^⊤ is a canonical vector in the vector space ℝ^N, in which the sth component is equal to one and others are zeros.
In a network with N nodes, to describe the relations between any two nodes, the Kronecker product ⊗ is adopted. Note that E(st)=e(s)⊗e^⊤(t) is a second-order canonical tensor in the space ℝ^N× N.
Given the intensity u_st
of the link between nodes s and t, the tensor U∈ℝ^N× N encodes all links  in the network, where

    U=∑_s,t=1^Nu_stE(st).

Next, we follow the tensorial notation in which a^i (i=1,2,…,N) represents the ith component of a contravariant vector (column vector) a∈ℝ^N
while a_j  (j=1,2,…,N) represents the jth component of the corresponding covariant vector (row vector) a^⊤∈ℝ^N.
We use (s) to indicate the sth vector, and (st) to denote the (st)-th tensor. Further,  i and j are used to represent the corresponding components of a certain vector or tensor.
The adjacency tensor U can be written as a 1-covariant and 1-contravariant tensor

    U_j^i=∑_s,t=1^Nu_ste^i(s)e_j(t)=∑_s,t=1^Nu_stE^i_j(st),

where E^i_j(st) indicates the component of E(st) in the ith row and jth column.


In multilayer networks, there may be various intra- and interlayer links. We further introduce the second-order adjacency tensor C^i_j(θϕ) (θ,ϕ=1,2,…,K)  to describe them.
To avoid confusion, Latin letters  indicate nodes in networks and Greek letters represent  layers in the multilayer networks, the same as described in the article.
Clearly, C^i_j(θθ) encodes the intralayer links in  layer θ, and C^i_j(θϕ) encodes the interlayer links from layer θ to layer  ϕ.
Further, u_st(θϕ) indicates the intensity from node s in  layer θ to node t in  layer ϕ.
Because of the varying influence between layers,  we use a canonical basis of layers.
The vector e^⊤(θ)∈ℝ^K(θ=1,2,…,K) is the θth covariant canonical vector  and e^α(θ) represents the αth component of e^⊤(θ), while e_β(θ) denotes the βth component of the corresponding contravariant canonical vector e_θ.
The set of all second-order tensors E^α_β(θϕ)=e^α(θ)e_β(ϕ) forms the canonical basis of the space ℝ^K× K.
The multilayer adjacency tensor M^iα_jβ can be represented by the tensor product of
C^i_j(θϕ) and E^α_β(θϕ):


    M^iα_jβ=   ∑_θ,ϕ=1^K C^i_j(θϕ)E^α_β(θϕ)
    
    =   ∑_θ,ϕ=1^K∑_s,t=1^Nu_st(θϕ)π^iα_jβ(stθϕ),

where the fourth-order tensor π^iα_jβ(stθϕ)=e^i(s)e_j(t)e^α(θ)e_β(ϕ) forms the canonical basis in the space ℝ^N× N× K× K.

Introducing  w_αβ to indicate the interlayer influence from  layer α to  layer β, the tensor W∈ℝ^K× K encodes the interlayer influence between any two layers in  multilayer networks, where W can be written as the 1-covariant and 1-contravariant tensor

    W^α_β=∑_θ,ϕ=1^K w_αβe^α(θ)e_β(ϕ)=∑_θ,ϕ=1^K w_αβE^α_β(θϕ).

Finally, we introduce the fourth-order tensor H, called the interaction tensor,  in terms of the  tensor product

    H^iα_jβ=   W^α_β· M^iα_jβ
    
    =   ∑_θ,ϕ=1^K∑_s,t=1^Nu_st(θϕ)w_αβπ^iα_jβ(stθϕ),

which encodes the interactions between any two nodes in the multilayer networks.

To simplify the notation, the Einstein summation convention is  adopted. Repeated indices in a term, where one index is a subscript and the other is a superscript, denote the summation  over all indices, which is also called  a contraction in tensorial operations. For example, we use the convention in the following equations:

    A_i^i   =Σ_i=1^NA_i^i,
    
    A_j^i B^j_i   =Σ_i=1^NΣ_j=1^NA_j^i B^j_i,
    
    A_jβ^iα B_iα   =Σ_i=1^NΣ_α=1^K A_jβ^iα B_iα.

This convention can be adopted for the product of any number of tensors with arbitrary order.



  §.§.§ II. Tensor decomposition

Tensors can be decomposed into lower-dimensional objects like matrices and supravectors<cit.>. We first flatten the adjacency tensor M into a matrix, i.e., we represent the fourth-order tensor M∈ℝ^N× N× K× K by a matrix M∈ℝ^NK× NK.
The overline here allows us to distinguish a tensor from its lower-dimensional form, such as a matrix.
Thus, the matrix M is a block matrix

    M= [ (@,20pt,20pt)c.c.c.cc.c.c.cM_11   M_12   ⋯   M_1K
    M_21   M_22   ⋯   M_2K
    ⋮   ⋮   ⋱   ⋮
    M_K1   M_K2   ⋯   M_KK ],

where the block M_αβ∈ℝ^N× N encodes the links from  layer α to layer β.

We can vectorize the second-order tensor Φ_iα∈ℝ^N× K into a supravector Φ∈ℝ^NK:

    Φ= [ (@,20pt,20pt)cc.c.c.cΦ_1
    Φ_2
    ⋮
    Φ_K ],

where Φ_α∈ℝ^N encodes the multicentrality of nodes in  layer α.

Denoting the matrix form of the influence tensor W by W∈ℝ^K× K, we have

    W= [ (@,20pt,20pt)c.c.c.cc.c.c.c
    w_11    w_12   ⋯    w_1K
    
    w_21    w_22   ⋯    w_2K
    ⋮   ⋮   ⋱   ⋮
    
    w_K1    w_K2   ⋯    w_KK ],

where element w_αβ quantifies the interlayer influence from layer α to  layer β.
Using these decompositions, the tensorial equation

    H_jβ^iαΦ_iα=λ_βΦ_jβ

becomes a  matrix equation


    [ (@,20pt,20pt)c.c.c.cc.c.c.c
    w_11·M_11    w_12·M_12   ⋯    w_1K·M_1K
    
    w_21·M_21    w_22·M_22   ⋯    w_2K·M_2K
    ⋮   ⋮   ⋱   ⋮
    
    w_K1·M_K1    w_K2·M_K2   ⋯    w_KK·M_KK ]^⊤·[ (@,20pt,20pt)cc.c.c.cΦ_1
    Φ_2
    ⋮
    Φ_K ]= [ (@,20pt,20pt)cc.c.c.cλ_1·Φ_1
    λ_2·Φ_2
    ⋮
    λ_K·Φ_K ].


We further define the vector Λ=[λ_1,λ_2,⋯,λ_K]^⊤∈ℝ^K. Then the Eq. <ref> can be written as

    (W⊙M)^⊤·Φ=Λ⊙Φ,

where ⊙ denotes the Khatri-Rao product.

In the calculation of the global multicentrality in multilayer networks, we take

    w_αβ=f(Φ_α)/f(Φ_β),

and we prove that there exists a unique normalized supravector Φ satisfying

    (W⊙M)^⊤·Φ=λ_1·Φ,

where λ_α=λ_1 (α=1,2,⋯,K).

In the calculation of the local multicentrality in multilayer networks, we take

    w_αβ=∑_i,j=1^NM^jβ_iαΦ_jβ/∑_i,j=1^NM^iα_jβΦ_iα,

and we also prove that there exists  a unique normalized supravector Φ satisfying

    (W⊙M)^⊤·Φ=Λ⊙Φ.


It is possible that different layers have different numbers of nodes, resulting in different tensor size. To maintain the unified tensor framework, we need to expand the tensor M and Φ to a specific size by adding zeros for nodes which are not in the current layer.
Thus, after tensor decomposition, many rows and the corresponding columns only have zero elements in the matrix M, and the corresponding multicentrality Φ_iα are also zeros. To reduce the computational complexity in practice, we can omit these rows and columns in the matrix M and the corresponding elements in the supravector Φ.



  §.§.§ III. PageRank centrality in a single-layer network

Consider a graph G=(V,E), where V={v_1,v_2,⋯,v_N} is the set of nodes and E is the set of  directed links. We have the adjacency matrix A=(a_ij), where

    a_ij=
    
    1,   if (i,j)∈ E
    
    0,   otherwise.

Further, we construct a stochastic matrix T=(t_ij), which is the transposition of the adjacency matrix A after normalization, i.e.,

    t_ij=
    1/k^out_j,   if (j,i)∈ E
    
    0,   otherwise

where k^out_j is the out-degree  of node j. The PageRank  centrality C∈ℝ^N<cit.> is the normalized leading right eigenvector of T, which corresponds to the leading eigenvalue λ_1=1, i.e., 

    T· C=1· C.


Next, we will briefly prove the existence and uniqueness of C.
We first assume that the graph G is strongly connected. 
Clearly, the summation of each column of T is equal to one, and T has a left eigenvector 1_1× N, whose components are all equal to one, corresponding to the case of eigenvalue λ=1.  Since ∑_i|t_ij|=1, all eigenvalues satisfy |λ_j|≤1 according to Gershgorin's circle theorem<cit.>.
Furthermore,  the matrix T is irreducible since the graph G is strongly connected. Hence, the leading eigenvalue λ_1 of T is unique  according to the Perron-Frobinus theorem<cit.>, leading to the fact that the normalized leading eigenvector C is unique.

When the graph G is not strongly connected, the following two steps can be employed to guarantee  the existence and uniqueness of the PageRank centrality. At the first step, find all dead ends (nodes with no outgoing links) in the graph G, and add outgoing links from them to all nodes.   Since we add outgoing links between a dead end and all nodes in the graph, such an operation will not change the relative centrality of other nodes<cit.>. At the second step, introduce a damping factor d ranging between zero and one<cit.>, and  the corresponding stochastic matrix T_d becomes

    T_d=d· T+1-d/N·1_N× N.

This makes  T_d irreducible, and thus the leading eigenvalue λ_1 of T is unique<cit.>.  

The PageRank centrality can be interpreted by two models: the flow model<cit.> and the random walk model<cit.>.
In the flow model,  the PageRank score is regarded as a kind of negotiable flow.
At each iteration, the current flow (score) of each node is distributed through its outgoing links on an equal footing to its neighboring nodes.
In the steady state, the flow that each node gains is its PageRank score. In the random walk model, we consider that a random walker  starts at a certain node and randomly chooses an outgoing link to walk iteratively.
In this context, the PageRank score of each node is the probability of the  walker resting on the node in the steady state.



  §.§.§ IV. PageRank multicentrality in multilayer networks

In multilayer networks, we define the stochastic tensor T_jβ^iα=M_jβ^iα/∑_j,βM_jβ^iα such that v^jβT_jβ^iα=v^iα, where v is the tensor with all components equal to one.
Further, we define the interaction tensor as H_jβ^iα=W^α_β· T_jβ^iα. It follows from Eq. <ref> that the leading eigentensor of H, denoted by Φ_iα, is  the PageRank multicentrality of node i in layer α, i.e.,

    H_jβ^iαΦ_iα=1·Φ_jβ.

Also, the influence tensor W can be a function of the tensor T and centrality Φ.
It is noteworthy that in PageRank multicentrality, the leading eigenvalue λ_1 of the interaction tensor is equal to one in accordance with PageRank centrality in a single-layer network.
As before, a damping factor d ∈ [0,1] should be introduced to guarantee the existence and uniqueness of global and local PageRank multicentrality when the multilayer networks are not strongly connected (see Text III of the SI above for further details).



Similarly, to characterize global PageRank multicentrality, we can choose the influence tensor as W^α_β=f(Φ_:α)/f(Φ_:β).
In particular, after a node in layer α equally distributes its multicentrality score to its neighbors, the scores propagating along intralayer links stay unchanged, while the scores propagating along interlayer links to nodes in layer β will be scaled by a factor f(Φ_:α)/f(Φ_:β).
In this case, we prove that the leading eigenvalue of H is equal to 1 (see Text V of the SI below for further details about the theoretic analysis).
For clarification, we provide  a toy example in Supplementary Fig. <ref>.


In local PageRank multicentrality, we specify the influence tensor as 
    W^α_β=∑_i,j=1^NT^jβ_iαΦ_jβ/∑_i,j=1^NT^iα_jβΦ_iα.

Likewise, we prove that the leading eigenvalue of H remains one (see Text VI of the SI below for the proof).
When a node in layer α propagates its multicentrality score to nodes in layer β, the interactions depend on the scores from nodes in layer β  to nodes in  layer α.
In other words, the scores leaving from layer α will be returned to layer α. We provide an illustration of local PageRank multicentrality in Supplementary Fig. <ref>.
In the random walk model, the interaction  tensor H^iα_jβ encodes the transition probability of node i in layer α leaving for any node j in layer β.
We can interpret the local multicentrality by a random walk process over a multilayer network as follows.
If a walker starts at a certain node in layer α, she will walk along the intralayer links normally according to the transition probability.
However, once the walker  leaves from  layer α to any other layer along the interlayer links, she will not rest on nodes in other layers until she goes back to layer α. The local PageRank multicentrality score of each node in layer α is the probability of such a walker resting on the node in the steady state.



  §.§.§ V. Theoretical demonstration of global multicentrality


Theorem 1: Given a strongly connected  multilayer network ℳ and a function f:ℝ^N_+→ℝ_+, we denote the eigenvector centrality in the corresponding single-layer network by Θ. Then, a unique global eigenvector multicentrality Φ exists if and only if, for ∀ Θ_:α, Θ_:β (α,β=1,2,⋯,K), there exists only one y_αβ∈ (0,+∞) such that f(y_αβ·Θ_:β)=1/y_αβ· f(Θ_:α).

Proof of Theorem 1.
(Sufficiency.) First, we will prove the existence of Φ.
Because the multilayer network is strongly connected,  there is a unique normalized leading eigentensor Θ (i.e., the eigenvector centrality in the corresponding single-layer network), satisfying

    M_jβ^iαΘ_iα=λ_1Θ_jβ,

which we have stated in Text III.
Then for Θ_:1∈ℝ^N and Θ_:γ∈ℝ^N (γ =1,2,⋯,K), there exists only one y_1γ∈ (0,+∞) such that f(y_1γ·Θ_:γ)=1/y_1γ· f(Θ_:1).
Let Φ_:γ=y_1γ·Θ_:γ (γ=1,2,⋯,K). We have

    ∑_α=1^K∑_i=1^NH^iα_jβΦ_iα=   ∑_α=1^K∑_i=1^N(w_αβ M^iα_jβ) · (y_1α·Θ_iα)
    
    =   y_1β∑_α=1^K∑_i=1^N M^iα_jβΘ_iα
    
    =   λ_1Φ_jβ.


Then, we will prove the uniqueness of Φ by contradiction. Suppose there exist two tensors Φ' and Φ” such that ∑_i,αΦ'_iα=∑_i,αΦ”_iα and

    H_jβ^iαΦ'_iα=λ_1Φ'_jβ,
    
    H_jβ^iαΦ”_iα=λ_1Φ”_jβ.

For Φ'_iα, we have

    ∑_i,αH_jβ^iαΦ'_iα   =λ_1·Φ'_jβ
    ∑_i,αf(Φ'_:α)/f(Φ'_:β)M_jβ^iαΦ'_iα   =λ_1·Φ'_jβ
    ∑_i,αM_jβ^iα· f(Φ'_:α)Φ'_iα   =λ_1· f(Φ'_:β)Φ'_jβ.

Similarly for Φ”_:α, we  have

    ∑_i,αM_jβ^iα· f(Φ”_:α)Φ”_iα   =λ_1· f(Φ”_:β)Φ”_jβ.

Due to the uniqueness of the leading eigentensor Θ of M, we have f(Φ'_:α)·Φ'_iα=f(Φ”_:α)·Φ”_iα=k·Θ_:α, and  Φ'=Φ”. This is in contradiction with our hypothesis. Thus, the solution Φ is unique.

(Necessity.) We proceed to prove the necessity of Theorem by contradiction.
Consider a strongly connected  multilayer network ℳ and a function f:ℝ^N_+→ℝ_+. Suppose that a unique global eigenvector multicentrality Φ exists. Further, we suppose that
∃ Θ_:μ,Θ_:ν, and either of the following two conditions is satisfied: 1) ∄ y_μν such that f(y_μν·Θ_:ν)=1/y_μν· f(Θ_:μ); 2)  ∃ y'_μν,y”_μν∈ (0,+∞) such that f(y'_μν·Θ_:ν)=1/y'_μν· f(Θ_μ) and f(y”_μν·Θ_ν)=1/y”_μν· f(Θ_μ).


By hypothesis, there is  a unique global eigenvector multicentrality Φ∈ℝ^N× K for the tensor M, such that

    H_jβ^iαΦ_iα=λ_1Φ_jβ,

i.e.,

    ∑_i,αH_jβ^iαΦ_iα   =λ_1·Φ_jβ
    ∑_i,αf(Φ_:α)/f(Φ_:β)M_jβ^iαΦ_iα   =λ_1·Φ_jβ
    ∑_i,αM_jβ^iα· f(Φ_:α)Φ_iα   =λ_1· f(Φ_:β)Φ_jβ.

Then we have


    f(Φ_:μ)·Φ_:μ=k·Θ_:μ
    
    f(Φ_:ν)·Φ_:ν=k·Θ_:ν.

Further, there exist y_μμ>0 and y_μν>0 such that y_μμ/y_μν is a constant, and the following two equations hold: 

    Φ_:μ=y_μμ·Θ_:μ
    Φ_:ν=y_μν·Θ_:ν.

Without loss of generality, we let y_μμ=1, i.e., Φ_:μ=Θ_:μ.
Thus, there exists only one y_μν>0 such that

    f(Φ_:μ)=f(Φ_:ν)· y_μν=k,

i.e.,

    f(y_μν·Θ_:ν)=1/y_μν· f(Θ_:μ),

which is in contradiction with the two conditions of the hypothesis.
▪


Note that Φ_:α represents the eigenvector multicentrality scores of nodes in layer α, Θ_:β represents the eigenvector centrality scores of nodes in layer β, and Φ_:β=y_αβ·Θ_:β represents the eigenvector multicentrality of nodes in layer β. Thus, y_αβ is a scale factor comparing the eigenvector multicentrality in a multilayer network with the eigenvector centrality in the corresponding single-layer network. If y_αβ<1, the eigenvector multicentrality score of nodes in layer β is relatively less important than the eigenvector centrality score in the corresponding single-layer network when taking layer α as a reference (Φ_:α=Θ_:α), and vice versa.
Hence, it is necessary to guarantee the existence and  uniqueness of y_αβ for any Θ_:α and Θ_:β as the theorem states. In fact, the condition in the theorem requires that the slope of the function f  has a lower bound. So, if the function f is a nonnegative monotonic increasing function such as the L^p-norm, the condition in the theorem is satisfied since the slope of  a  monotonic increasing function is always positive.


In addition, from the proof above, we can see that λ_1 is the leading eigenvalue of the tensor M, satisfying

    M_jβ^iαΘ_iα=λ_1Θ_jβ.

Also, λ_1 is the leading eigenvalue of the tensor H, satisfying

    H_jβ^iαΦ_iα=λ_1Φ_jβ.


In the global PageRank multicentrality, the leading eigenvalue of the tensor M is equal to one. Clearly, the leading eigenvalue of the tensor H is also equal to one.



  §.§.§ VI. Theoretical demonstration of local multicentrality


Theorem 2: Given a strongly connected multilayer network ℳ and influence matrix W described by

    w_αβ=∑_i,j=1^NM^jβ_iαΦ_jβ/∑_i,j=1^NM^iα_jβΦ_iα,

the  local eigenvector multicentrality Φ has a unique normalized solution to the equation

    H_jβ^iαΦ_iα=λ_βΦ_jβ.


Proof of Theorem 2.
We  flatten the interaction tensor H into a matrix H∈ℝ^NK× NK, and vectorize the second-order tensor Φ∈ℝ^N× K into a vector Φ∈ℝ^NK.
We will show that  H·Φ is a contraction mapping of Φ. To this end, we first suppose that the adjacency tensor M satisfies v^jβM_jβ^iα=v^iα, where v is the tensor with all components equal to one, i.e., the case of local PageRank multicentrality.
Consider two different vectors Φ'_γ and Φ”_γ. For any given Φ_η (η≠γ), we suppose

    Φ'(γ)= [ (@,20pt,20pt)cc.c.c.c.c.cΦ_1
    Φ_2
    ⋮
    Φ'_γ
    ⋮
    Φ_M ]∈ℝ^N,  Φ”(γ)= [ (@,20pt,20pt)cc.c.c.c.c.cΦ_1
    Φ_2
    ⋮
    Φ”_γ
    ⋮
    Φ_M ]∈ℝ^N.

Then



    |H·Φ'(γ)-H·Φ”(γ)|
    
    =   |M_γγΦ'_γ+∑_η≠γM_ηγ·Φ'_η_1/M_γη·Φ'_γ_1·M_γη·Φ'_γ-M_γγΦ”_γ-∑_η≠γM_ηγ·Φ”_η_1/M_γη·Φ”_γ_1·M_γη·Φ”_γ|
    ⩽   p_γ· |Φ'_γ-Φ”_γ|,

where p_γ is given by

    p_γ=max_i,k∑_jM_γγ(i,j)-M_γγ(k,j)_1<1.

Thus, for any positive vectors Φ' and Φ”,

    |H·Φ'-H·Φ”|⩽max_γ{p_γ}· |Φ'-Φ”|.

This implies that the function H·Φ is a contraction mapping of Φ. According to the Banach fixed-point theorem<cit.>, H admits a unique fixed-point Φ, satisfying

    H·Φ=Φ,

from which we find that the eigenvalue of H_jβ^iα remains one in the local PageRank multicentrality analysis.
For the general case, we can construct an analogous contraction mapping.
▪



  §.§.§ VII. Epidemic spreading in  multilayer networks


We first introduce epidemic spreading process in a single-layer network. There are many models describing epidemic  spreading process, such as the susceptible-infected (SI) model, the susceptible-infected-susceptible (SIS) model, the susceptible-infected-recovered (SIR) model and so on<cit.>. Here we take the SI model as an illustration.

In an undirected single-layer  social network, nodes denote individuals and links indicate their social relations. We study the spreading of a disease in such a social network.
Each node of the network has only  two states, “susceptible” (S) or “infected” (I). A susceptible node is a temporarily healthy node, which can be infected by infected nodes, while an infected node will stay in the infected state.
At each step, an infected node will infect the susceptible nodes connected to it with probability μ, which is called the transition rate. Clearly, if a susceptible node has more infected neighbors, it will have a higher probability of being infected.
Since the infected nodes will not  recover, all nodes will eventually  be infected  if the network is strongly connected and μ>0.

In multilayer networks, nodes lie in different layers. The transition rate can be different along intralayer links and interlayer links. Here we consider a dataset  from a large European research institute<cit.> (data provided by the Stanford Network Analysis Project, <http://snap.stanford.edu/index.html>). There are 1,005 individuals, each belonging to exactly one of 42 departments at the research institute. A multilayer network is built with  K=42 layers indicating departments, and N=1,005 nodes representing individuals. There is an undirected link between two nodes if either one of them has sent at least one email to the other, and the connections are encoded by the adjacency tensor M ∈ℝ^N× N× K× K.
We then consider an epidemic spreading process in this multilayer network.
Considering a certain node as the infection source, we set the intralayer transition rate to 0.05 and the interlayer transition rate to 0.03 at each  step<cit.>.  Since we already have prior knowledge about interlayer influence, the influence tensor W ∈ℝ^K× K is given as follows:

    W^α_β=
    
    0.05,   α=β
    
    0.03,   α≠β.


Then, we obtain the interaction tensor H ∈ℝ^N× N× K× K, and the multicentrality tensor Φ∈ℝ^N× K by solving the tensorial equation

    H_jβ^iαΦ_iα=λ_1Φ_jβ.



We simulate the spreading process in this multilayer network with only one infected node in the initial state. At each  step, the infected node/nodes  will infect its/their susceptible neighboring nodes at the transition rate, and the rates are different between intralayer links and interlayer links as we mentioned above.
All nodes will be infected as long as there is a path between any two nodes in the multilayer network. We focus on the spreading coverage, i.e., the percentage of the infected nodes after each step.
The results are shown in Supplementary Fig. <ref>, indicating that the nodes with higher eigenvector multicentrality play a more important role in the process of epidemic spreading in a multilayer network.



  §.§.§ VIII. Theoretical demonstration of the numerical solution


Here we prove that the iteration scheme in the Materials and Methods section,

    Φ^(k+1)=Φ^(k)+D^(k)[(W⊙M)^⊤(k)⊙Ω^(k)-E_N]Φ^(k),

converges to the solution to the matrix equation

    (W⊙M)^⊤·Φ=Λ⊙Φ.


Proof. We first consider the global multicentrality. The iteration scheme Eq. <ref> becomes

    Φ^(k+1)=Φ^(k)+D^(k)[(W⊙M)^⊤(k)/λ_1^(k)-E_N]Φ^(k),

and the matrix Eq. <ref> becomes

    (W⊙M)^⊤·Φ=λ_1⊙Φ.



In Text V of the SI above, we have proved that there exists a unique solution to Eq. <ref>.
We denote the unique solution by Φ^* and the corresponding eigenvalue by λ_1^*.
The convergence of the iteration scheme Eq. <ref> has been demonstrated in the Materials and Methods section. We denote the limit of Φ^(k) by Φ^(∞). The primary task is to prove that  Φ^(∞) is equal to Φ^*.
Clearly,

    Φ^(∞)   =Φ^(∞)+D^(∞)[(W⊙M)^⊤(∞)/λ_1^(∞)-E_N]Φ^(∞)
    
    0   =D^(∞)[(W⊙M)^⊤(∞)/λ_1^(∞)-E_N]Φ^(∞).

As the matrix D^(∞) is a diagonal matrix with non-zero diagonal elements, it is invertible. Thus

    [(W⊙M)^⊤(∞)/λ_1^(∞)-E_N]Φ^(∞)   =0
    
    (W⊙M)^⊤(∞)/λ_1^(∞)·Φ^(∞)   =Φ^(∞).

Therefore, we have

    (W⊙M)^⊤(∞)·Φ^(∞)=λ_1^(∞)·Φ^(∞).

Note that Φ^(∞) satisfies Eq. <ref> and has a unique solution Φ^*. Thus,  we have Φ^(∞)=Φ^* and λ_1^(∞)=λ_1^*. Similarly, we can prove the convergence of the iteration scheme for local multicentrality. ▪



 §.§ Supplementary Figures

































 §.§ Supplementary Tables




































































§ ACKNOWLEDGEMENTS.

We thank Manlio De Domenico for sharing the processed dataset of the European airline network and the details for calculating the versatility. This work was  supported by the National Natural Science Foundation of China (Grant No. 61731004) and the Zhejiang Natural Science Foundation (Grant No. LR16F020001).  This work was supported in part by the U.S. Army Research Office under Grant W911NF-16-1-0448 and DTRA under Grant HDTRA1-13-1-0029.



§ AUTHOR CONTRIBUTIONS.

S.H., J.C., J.Z., and H.V.P. conceived and designed the project.
M.W. conducted all the theoretical analysis and empirical calculations.
All authors analyzed the results.
M.W., S.H., J.Z., and J.C. wrote the manuscript.
Y.-Y.L., J.Z., H.V.P., and Y.S. edited the manuscript.





 §.§ Competing Interests.

The authors declare no competing financial interests.





§ REFERENCES
















