Tunable Eigenvector-Based Centralities for Multiplex and Temporal Networks

TL;DR

This paper introduces a unified linear-algebraic framework for eigenvector-based centralities in multiplex and temporal networks, accounting for complex interlayer couplings and their influence on node importance measures.

Contribution

It generalizes existing centrality measures to multilayer networks with asymmetric interlayer couplings, extending prior work to include a broader class of network interactions.

Findings

Framework unifies centrality analysis for multiplex and temporal networks.

Interlayer coupling topology significantly affects node centralities.

Singular perturbation theory describes limits of weak and strong interlayer coupling.

Abstract

Characterizing the importances (i.e., centralities) of nodes in social, biological, and technological networks is a core topic in both network science and data science. We present a linear-algebraic framework that generalizes eigenvector-based centralities, including PageRank and hub/authority scores, to provide a common framework for two popular classes of multilayer networks: multiplex networks (which have layers that encode different types of relationships) and temporal networks (in which the relationships change over time). Our approach involves the study of joint, marginal, and conditional "supracentralities" that one can calculate from the dominant eigenvector of a supracentrality matrix [Taylor et al., 2017], which couples centrality matrices that are associated with individual network layers. We extend this prior work (which was restricted to temporal networks with layers that are coupled by adjacent-in-time coupling) by allowing the layers to be coupled through a (possibly asymmetric) interlayer-adjacency matrix $\tilde{{\bf A}}$, where the entry $\tilde{A}_{tt'} \geq 0$ encodes the coupling between layers $t$ and $t'$. Our framework provides a unifying foundation for centrality analysis of multiplex and temporal networks; it also illustrates a complicated dependency of the supracentralities on the topology and weights of interlayer coupling. By scaling $\tilde{{\bf A}}$ by an interlayer-coupling strength $\omega\ge0$ and developing a singular perturbation theory for the limits of weak ($\omega\to0^+$) and strong coupling ($\omega\to\infty$), we also reveal an interesting dependence of supracentralities on the dominant left and right eigenvectors of $\tilde{{\bf A}}$.

Figures

Tables

Table 1. Table 1: Summary of our mathematical notation for objects of different dimensions.

Typeface	Class	Dimension
$𝕄$	matrix	$NT\times NT$
$𝐌$	matrix	$N\times N$
$𝑴$	matrix	$T\times T$
$𝕧$	vector	$NT\times 1$
$𝐯$	vector	$N\times 1$
$𝒗$	vector	$T\times 1$
$M_{ij}$	scalar	1
$v_i$	scalar	1

Table 2. Table 2: A list of European airline companies, which we represent as layers in a multiplex network. For each layer, we report the number M l subscript 𝑀 𝑙 M_{l} of undirected edges and the spectral radius λ 1 ( t ) superscript subscript 𝜆 1 𝑡 \lambda_{1}^{(t)} of its associated adjacency matrix 𝐀 ( t ) superscript 𝐀 𝑡 {\bf A}^{(t)} . We have chosen the ordering to match the one in [ 13 ] .

Layer $t$ Airline Name $M_t$ ${\lambda }_1^{\left(t\right)}$ 1 Lufthansa 244 14.5 2 Ryanair 601 19.3 3 easyJet 307 14.0 4 British Airways 66 6.6 5 Turkish Airlines 118 9.9 6 Air Berlin 184 11.3 7 Air France 69 7.2 8 Scandinavian Air. 110 8.9 9 KLM 62 7.9 10 Alitalia 93 8.8 11 Swiss Int. Air Lines 60 7.3 12 Iberia 35 5.8 13 Norwegian Air Shu. 67 8.1 14 Austrian Airlines 74 8.1 15 Flybe 99 8.5 16 Wizz Air 92 6.5 17 TAP Portugal 53 7.0 18 Brussels Airlines 43 6.6 19 Finnair 42 6.4

Layer $t$ Airline Name $M_t$ ${\lambda }_1^{\left(t\right)}$ 20 LOT Polish Air. 55 6.8 21 Vueling Airlines 63 6.8 22 Air Nostrum 69 6.4 23 Air Lingus 108 6.7 24 Germanwings 67 7.4 25 Pegasus Airlines 58 6.7 26 NetJets 180 8.2 27 Transavia Holland 57 6.0 28 Niki 37 4.7 29 SunExpress 67 7.8 30 Aegean Airlines 53 6.5 31 Czech Airlines 41 6.4 32 European Air Trans. 73 6.8 33 Malev Hungarian Air. 34 5.8 34 Air Baltic 45 6.4 35 Wideroe 40 5.6 36 TNT Airways 61 6.2 37 Olympic Air 43 6.2

Table 3. Table 3: European airports with the largest MNCs for interlayer-coupling strengths ω 𝜔 \omega in the regimes of weak ( ω = 0.01 𝜔 0.01 \omega=0.01 ), intermediate ( ω = 1 𝜔 1 \omega=1 ), and strong ( ω = 100 𝜔 100 \omega=100 ) coupling. We give the results of our computations of eigenvector supracentralities with all-to-all, uniform, and diagonal coupling between layers. We identify each airport by its International Civil Aviation Organization (ICAO) code.

Rank 1 2 3 4 5 6 7 8 9 10

$\omega =0.01$ Airport MNC EGSS 0.329 EIDW 0.286 LIME 0.254 EBCI 0.201 LEMD 0.193 LEAL 0.190 EDFH 0.189 LIRA 0.184 LEGE 0.176 LEPA 0.166

$\omega =1$ Airport MNC LEMD 1.379 EHAM 1.296 LEBL 1.257 EDDM 1.171 LIRF 1.150 EDDF 1.121 EDDL 1.105 LFPG 1.091 LOWW 1.066 LIMC 0.968

$\omega =100$ Airport MNC EHAM 1.406 LEMD 1.400 LIRF 1.206 LOWW 1.198 LEBL 1.193 EDDM 1.160 LFPG 1.157 EDDF 1.134 EDDL 1.128 LSZH 1.017

Equations

$$\mathbf{Cu}={\lambda }_{\max }\mathbf{u},$$

$$\mathbb{C}(\omega )$$

$$\mathbb{C}(\omega )\text{mathbbm}v(\omega )={\lambda }_{\max }(\omega )\text{mathbbm}v(\omega ),$$

$$W_{it}(\omega )=\text{mathbbm}{v}_{N(t-1)+i}(\omega ).$$

$$x_t(\omega )=\sum _i{W}_{it}(\omega ),{\hat{x}}_i(\omega )=\sum _t{W}_{it}(\omega ).$$

$$Z_{it}(\omega )=W_{it}(\omega )/x_t(\omega ),{\hat{Z}}_{it}(\omega )=W_{it}(\omega )/{\hat{x}}_i(\omega ),$$

$$\mathbb{C}(0)\left(\sum _{(i,t)\in \mathcal{P}}{\alpha }_{i,t}\text{mathbbm}{v}^{(i,t)}\right)=\sum _{(i,t)\in \mathcal{P}}{\alpha }_{i,t}\mathbb{C}(0)\text{mathbbm}{v}^{(i,t)}$$

$$=\lambda \left(\sum _{(i,t)\in \mathcal{P}}{\alpha }_{i,t}\text{mathbbm}{v}^{(i,t)}\right).$$

$$\mathbbm{v}^{(1)}(\omega)\to\sum_{t\in{{\color[rgb]{0,0,0}\mathcal{T}}}}\alpha_{t}\mathbbm{v}^{(1,t)}\,,~{}~{}\mathbbm{u}^{(1)}(\omega)\to\sum_{t\in{{\color[rgb]{0,0,0}\mathcal{T}}}}\beta_{t}\mathbbm{u}^{(1,t)}\,,$$

$$\bm{X}\bm{\alpha}=\lambda_{1}\bm{\alpha}\,,~{}~{}{\bm{X}}^{{\color[rgb]{0,0,0}*}}\bm{\beta}=\lambda_{1}\bm{\beta}\,,$$

$$X_{t{t}^{\prime }}$$

$$\displaystyle\mathbbm{v}(\omega)\to\mathbbm{v}^{(1,{{\color[rgb]{0,0,0}\tau}})}\,,~{}~{}\mathbbm{u}(\omega)\to\mathbbm{u}^{(1,{{\color[rgb]{0,0,0}\tau}})}\,.$$

$$\tilde{\mathbb{C}}(ϵ)=ϵ\mathbb{C}({ϵ}^{-1})=ϵ\hat{\mathbb{C}}+\hat{\mathbb{A}},$$

$$[\mathbb{P}]_{kl}=\left\{\begin{array}[]{rl}1\,,&l=\lceil k/N\rceil+T\,[(k-1)\bmod N]\\ 0\,,&\textrm{otherwise}\,.\end{array}\right.$$

$${\widetilde{\text{mathbbm}v}}^{(t)}=\sum _j{\tilde{\alpha }}_{tj}\mathbb{P}{\widetilde{\text{mathbbm}v}}^{(t,j)},{\widetilde{\text{mathbbm}u}}^{(t)}=\sum _j{\tilde{\beta }}_{tj}\mathbb{P}{\widetilde{\text{mathbbm}u}}^{(t,j)},$$

$${\tilde{\lambda }}_{\max }(ϵ)\to {\tilde{\mu }}_1,{\widetilde{\text{mathbbm}v}}^{(1)}(ϵ)\to \sum _j{\tilde{\alpha }}_j\mathbb{P}{\widetilde{\text{mathbbm}v}}^{(1,j)},{\widetilde{\text{mathbbm}u}}^{(1)}(ϵ)\to \sum _j{\tilde{\beta }}_j\mathbb{P}{\widetilde{\text{mathbbm}u}}^{(1,j)},$$

$$\tilde{\mathbf{X}}\tilde{\bm{\alpha}}=\tilde{\mu}_{1}\tilde{\bm{\alpha}}\,,~{}~{}\tilde{\mathbf{X}}^{{\color[rgb]{0,0,0}*}}\tilde{\bm{\beta}}=\tilde{\mu}_{1}\tilde{\bm{\beta}}\,,$$

$${\tilde{X}}_{ij}$$

$$\tilde{\mathbf{X}}=\sum _t{\mathbf{C}}^{(t)}\frac{{\sin }^2\left(\frac{\pi t}{T+1}\right)}{{\sum }_{t=1}^T{\sin }^2\left(\frac{\pi t}{T+1}\right)}.$$

$$\tilde{A}_{tt^{\prime}}=\left\{\begin{array}[]{rl}1\,,&t^{\prime}-t=1\\ \gamma\,,&\text{otherwise}\,.\\ \end{array}\right.$$

$$\{(k_0,t_0),(k_1,t_1),(k_2,t_2),\dots ,(k_{l-1},t_{l-1}),(k_l,t_l)\},$$

$$\displaystyle[{{\color[rgb]{0,0,0}\mathbb{C}(0)}}\mathbbm{v}^{(i,t)}]_{{\color[rgb]{0,0,0}p}}=\sum_{(k,t^{\prime})}C_{jk}^{(t)}v_{k}^{(i,t^{\prime})}\delta_{tt^{\prime}}=\mu_{i}^{(t)}v_{j}^{(i,t^{\prime})}\delta_{tt^{\prime}}=\mu_{i}^{(t)}[\mathbbm{v}^{(i,t)}]_{{\color[rgb]{0,0,0}p}}\,,$$

$${\lambda }_{\text{max}}(\omega )$$

$$\text{mathbbm}{v}^{(1)}(\omega )$$

$$\text{mathbbm}{u}^{(1)}(\omega )$$

$$\left({\lambda }_0\mathbb{I}-\hat{\mathbb{C}}\right)\text{mathbbm}{v}_0$$

$$\left({\lambda }_0\mathbb{I}-\hat{\mathbb{C}}\right)\text{mathbbm}{v}_1$$

$$\text{mathbbm}{v}_0=\sum _t{\alpha }_t\text{mathbbm}{v}^{(1,t)},\text{mathbbm}{u}_0=\sum _t{\beta }_t\text{mathbbm}{u}^{(1,t)},$$

$$\displaystyle[\mathbbm{u}^{(1,t)}]^{{\color[rgb]{0,0,0}*}}\hat{\mathbb{A}}\mathbbm{v}_{0}$$

$$\displaystyle\sum_{t^{\prime}}\alpha_{t^{\prime}}[\mathbbm{u}^{(1,t)}]^{{\color[rgb]{0,0,0}*}}\hat{\mathbb{A}}\mathbbm{v}^{(1,t^{\prime})}$$

Full text

{AMS}

91D30, 05C81, 94C15, 05C82, 15A18

1 Introduction

Quantifying the importance of entities in a network is an essential feature of many search engines on the World Wide Web [11, 80, 62, 33], ranking algorithms for sports teams and athletes [12, 96, 15], targeted social-network marketing schemes [53], investigations of fragility in infrastructures [43, 40], quantitative analysis of the impact of research papers and scientists [31], examinations of the influence of judicial and legislative documents [63, 32], identification of novel drug targets in biological systems [49], and many other applications. In the most common (and simplest) type of network, called a “graph” or a “monolayer network”, a node represents an entity (e.g., a web page, a person, a document, or a protein), and an edge encodes a relationship between a pair of entities. Centrality analysis, in which one seeks to quantify the importances of nodes and/or edges (and, more generally, of other subgraphs as well), has been developed intensively across numerous domains, including sociology, mathematics, computer science, and physics [19, 78, 33, 62].

Researchers have developed increasingly comprehensive network representations and analyses [85, 61] to help with data integration and a variety of applications. A prominent example is the generalization of graphs to multilayer networks [55, 8, 18, 84], and there have been many efforts to extend centrality measures to multiplex and temporal networks [42, 98, 99, 104, 2, 82, 54, 121, 64, 73, 39, 27, 38, 103, 93, 95, 123, 110]. Multilayer network centralities have been used in the study of diverse applications, including social networks [68, 67, 42, 17, 14], transportation systems [22, 105, 48, 116], economic systems [23, 24, 5], neural systems [6, 20, 48, 124], and signal processing of geological time series [66]. Moreover, many techniques in centrality analysis are closely connected to the study of various dynamical processes (including in multilayer networks), such as random walks [33, 72, 79, 35, 42, 25], information spreading [17, 91], and congestion [22].

We consider two types of multilayer networks (see Fig. 1): (1) multiplex networks, in which layers represent different types of relationships; and (2) temporal networks, in which layers represent different time instances or time periods. We extend the mathematical framework of supracentrality matrices, which we developed recently [110] to generalize eigenvector-based centralities (e.g., PageRank [11, 80, 33], eigenvector centrality [10], and hub and authority scores [57]) to multilayer representations of discrete-time temporal networks. Our approach involves coupling centrality matrices that are associated with individual layers into a larger supracentrality matrix and studying its dominant eigenvector111Technically, we study the eigenvector that is associated with the largest positive eigenvalue $\lambda_{\textrm{max}}$ of an irreducible nonnegative matrix. Because the other eigenvalues have magnitudes that are less than or equal to $\lambda_{\textrm{max}}$, we refer to this eigenvalue and its eigenvectors as “dominant”. to obtain joint, marginal, and conditional centralities (see Sec. 3.2) to quantify the importances of nodes, layers, and node-layer pairs. In this article, we generalize the supracentrality framework of [110] to multiplex networks, which integrate data sets that encode different types of relationships by coupling them as layers of a single multilayer network.

Generalizing centrality measures to multiplex networks and temporal networks are active areas of research [55, 45, 44, 8, 18] (see our discussion in Sec. 2.3), and our supracentrality framework is relevant for such efforts.222In principle, one can also use supracentrality matrices of higher dimensionality to study networks that are both multiplex and temporal, but we do not study any such examples in this paper. Our original formulation of supracentrality in [110] focused on temporal networks (see [111] for our more recent work), and it assumed a specific type of multilayer representation with adjacent-in-time coupling. We now extend supracentrality matrices to a broader class of multilayer networks by coupling layers via an interlayer-adjacency matrix $\tilde{\bm{A}}$, where $\tilde{{A}}_{tt^{\prime}}\geq 0$ encodes the (possibly asymmetric) coupling between layers $t$ and $t^{\prime}$. We assume “diagonal” interlayer coupling (see Sec. 2.2 and [55]), as we only connect instantiations of the same entity (i.e., node) across different layers. We also assume that the interlayer coupling is “uniform”, so all edges between layers $t$ and $t^{\prime}$ have the same weight $\omega\tilde{A}_{tt^{\prime}}\geq 0$. Multilayer networks with both diagonal and uniform interlayer coupling are said to be “layer-coupled” [55]. The value of $\omega$ determines how strongly the layers influence each other. We will show that choices for $\tilde{\bm{A}}$ and $\omega$ significantly affect supracentralities and are useful “tuning knobs” to consider when calculating and interpreting supracentralities.

To gain insight into the effects of $\tilde{\bm{A}}$ and $\omega$, we use singular perturbation theory to analyze the dominant eigenspace of supracentrality matrices in the limits of weak ($\omega\to 0^{+}$) and strong ($\omega\to\infty$) coupling. We show that these limits yield layer decoupling and a type of layer aggregation (which is a form of data fusion), respectively. There are many scenarios in which one couples matrices into a larger “supramatrix” — including the detection of multilayer community structure using a supramodularity matrix [74, 119] and the study of random walks and diffusion on multilayer networks via supra-Laplacian matrices [35, 89] — and our perturbative approach reveals insights about the utility of matrix coupling as a general technique for multimodal data integration. Specifically, the singular perturbation theory of Sec. 4 makes no explicit assumption that the block-diagonal matrices are centrality matrices, so our results also characterize the dominant eigenspaces of layer-coupled matrices in other applications, including ones that are unrelated to networks.

Our results in Sec. 4 characterize the decoupling and aggregation limits of supracentrality matrices. We illustrate that the limiting dominant eigenspace of a supracentrality matrix depends on a complicated interplay between many factors, including (i) the dominant eigenvectors of the centrality matrix of each layer; (ii) the dominant eigenvectors of the interlayer-adjacency matrix; and (iii) the spectral radii of the layers’ centrality matrices. For the $\omega\to\infty$ limit, the dominant eigenspace of a supracentrality matrix depends on a weighted average of the layers’ centrality matrices, with weights that are related to the dominant eigenspace of $\tilde{\bm{A}}$. A key factor for the $\omega\to 0^{+}$ limit is whether the layers’ individual centrality matrices have identical or different spectral radii. In the latter scenario, we identify and characterize an eigenvector-localization phenomenon in which one or more layers dominate the decoupling limit. Our layer-aggregation and decoupling limits are reminiscent of prior research on supra-Laplacian matrices [101, 35], but our coupling matrices and qualitative results both differ from such prior work.

We illustrate our framework with applications to two empirical, multimodal network data sets. First, we study the importances of European airports in a multiplex network in which layers represent different airlines [13]. We find, for example, that supracentralities in the weak-coupling limit are dominated by the Ryanair layer; among all layers, this one has the most edges, and its associated adjacency matrix has the largest spectral radius. For intermediate coupling strengths, we observe a centrality “boost” (i.e., an increase relative to the centralities of other nodes) for airports that are central both with respect to the Ryanair layer and with respect to a network that is associated with an aggregation of the network layers (specifically, the one from summing the layers’ adjacency matrices). We study these phenomena by comparing marginal node centralities to the nodes’ intralayer degrees and total degrees (which quantify, respectively, the number of edges of a node in each layer and its total number of edges across all layers). Our second focal example, which we construct using data from the Mathematics Genealogy Project [88, 107], is a temporal network that encodes the graduation and hiring of mathematicians at mathematical-science Ph.D. programs in the United States [75, 110]. Extending [110], we explore the effects of causality by implementing time-directed interlayer coupling along with layer teleportation (see Fig. 1(b) and our discussion in Sec. 5.2), which we define analogously to node teleportation in PageRank [33]. (See our recent book chapter [111] for further exploration of layer teleportation.) Similar to previous findings for causality-respecting centralities [29, 37], our approach boosts the centralities of node-layer pairs whose edges occur earlier in time (allowing them to causally influence more nodes). In the present paper, this phenomenon manifests as a boost in marginal layer centrality for older time layers. Our numerical experiments highlight the importance of exploring a diverse set of interlayer-coupling architectures $\tilde{\bm{A}}$ and strengths $\omega$ to identify application-appropriate parameter choices.

Our paper proceeds as follows. In Sec. 2, we present background information on eigenvector-based centralities, multiplex and temporal networks, and generalizing centralities for such networks. In Sec. 3, we present our supracentrality framework. In Sec. 4, we analyze the weak-coupling and strong-coupling limits. In Sec. 5, we study the two empirical data sets. We conclude in Sec. 6 and give the proofs for our main mathematical results in appendices. We present further numerical investigations in Supplementary Materials.

2 Background Information

We now discuss eigenvector-based centralities in Sec. 2.1, multiplex and temporal networks in Sec. 2.2, and extensions of eigenvector-based centralities for multiplex and temporal networks in Sec. 2.3.

2.1 Eigenvector-Based Centrality Measures

We start with a definition.

Definition 2.1** (Monolayer Network).**

*Let $\mathcal{G}(\mathcal{V},\mathcal{E})$ be a monolayer (i.e., single-layer) network with nodes $\mathcal{V}=\{1,2,\dots,N\}$ and a set $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}\times\mathbb{R}_{+}$ of positively-weighted edges, where $(i,j,w_{ij})\in\mathcal{E}$ if and only if there exists an edge from $i$ to $j$ with weight $w_{ij}$. We also encode this network (which is a weighted graph) by an $N\times N$ adjacency matrix $\mathbf{A}$ with entries ${A}_{ij}=w_{ij}$ if $(i,j,w_{ij})\in\mathcal{E}$ and ${A}_{ij}=0$ otherwise. *

One of the most popular approaches for quantifying the importances of nodes $\mathcal{V}$ in a network is to calculate the dominant eigenvector of a network-related matrix and interpret the eigenvector’s entries as a proxy for importance [78, 33].

Definition 2.2** (Eigenvector-Based Centrality).**

Let $\mathbf{C}=C(\mathbf{A})$ be a “centrality matrix” that we obtain via some function $C:\mathbb{R}^{N\times N}\mapsto\mathbb{R}^{N\times N}$ of the adjacency matrix $\mathbf{A}$ for a network $\mathcal{G}(\mathcal{V},\mathcal{E})$. Consider the right eigenvector $\mathbf{u}$, which is the solution to

[TABLE]

where $\lambda_{\rm{max}}\in\mathbb{R}$ is the largest positive eigenvalue of $\mathbf{C}$. Each entry $u_{i}$ specifies the eigenvector-based centrality that is associated with the function $C$ for node $i\in\mathcal{V}$. We refer to this eigenvalue and its associated eigenvectors as “dominant”.

Different choices for $C(\mathbf{A})$ yield different notions of centrality, and some are more useful than others. The following are among the most popular eigenvector-based centralities.

Definition 2.3** (Eigenvector Centrality [10]).**

*With the choice $\mathbf{C}^{(EC)}=\mathbf{A}$, Eq. (1) yields eigenvector centralities $\{u_{i}^{(EC)}\}$. *

Definition 2.4** (Hub and Authority Scores [57]).**

*With the choices \mathbf{C}^{(HS)}=\mathbf{A}\mathbf{A}^{{\color[rgb]{0,0,0}*}} and \mathbf{C}^{(AS)}=\mathbf{A}^{{\color[rgb]{0,0,0}*}}\mathbf{A}, Eq. (1) yields hub scores $\{u_{i}^{(HS)}\}$ and authority scores $\{u_{i}^{(AS)}\}$, respectively. The symbol $*$ denotes the transpose operator333In this paper, we define $T$ to be the number of layers, so we use the notation $*$ for the transpose operator to avoid confusion. The symbol $*$ is often used to indicate the conjugate transpose operator. This is also a correct interpretation in the present paper, because the matrices that we study have real-valued entries.. *

Remark 2.5**.**

*Hub scores and authority scores are, respectively, the left and right dominant singular vectors of $\mathbf{A}$. *

Definition 2.6** (PageRank [80, 33]).**

*Consider the choice \mathbf{C}^{(PR)}={\sigma}{(\mathbf{D}^{-1}\mathbf{A})^{{\color[rgb]{0,0,0}*}}}+(1-{\sigma})N^{-1}\mathbf{1}\mathbf{1}^{{\color[rgb]{0,0,0}*}}, where $\mathbf{D}$ is the diagonal matrix with entries $D_{ii}=\sum_{j}A_{ij}$. The quantity ${\sigma}\in[0,1)$ is a “node” teleportation parameter (we will assume that ${\sigma}=0.85$ in this paper), and $\mathbf{1}$ is a length-$N$ vector of ones (such that \mathbf{1}\mathbf{1}^{{\color[rgb]{0,0,0}*}} is an $N\times N$ matrix of ones). Using this choice in Eq. (1), we obtain PageRank centralities $\{u_{i}^{(PR)}\}$. *

Remark 2.7**.**

*Our definition of PageRank assumes that each edge has at least one out-edge. If there exist “dangling nodes” (which lack out-edges), one can add self-edges (i.e., either just to such nodes or to all the nodes). *

Remark 2.8**.**

*It is also common to compute PageRank centralities from a left eigenvector [33]. In the present paper, we use a right-eigenvector formulation to be consistent with the other eigenvector-based centralities. One can recover the other formulation by taking the transpose of the matrix $\mathbf{C}$ in Eq. (1). *

Remark 2.9**.**

*There are other possible teleportation strategies for PageRank, such as ones with a local bias or emphasis on other features [33]. In such cases, one replaces the matrix \mathbf{1}\mathbf{1}^{{\color[rgb]{0,0,0}*}} by \mathbf{u}\mathbf{1}^{{\color[rgb]{0,0,0}*}}, where the vector $\mathbf{u}$ encodes the biases. *

In most applications, it is important to choose the function $C$ to ensure that centralities are unique and strictly positive. It is common to use the following two theorems to guarantee these important features.

Theorem 2.10** (Perron–Frobenius Theorem for Nonnegative Matrices [4]).**

*Let $\mathbf{C}\in\mathbb{R}^{N\times N}$ be an irreducible square matrix with nonnegative entries. It follows that $\mathbf{C}$ has a simple largest positive eigenvalue $\lambda_{\rm{max}}$ and that its right and left eigenvectors are positive and unique. Moreover, if $\mathbf{C}$ is aperiodic, then $\lambda_{\rm{max}}>|\lambda_{i}|$ for any $\lambda_{i}\not=\lambda_{\rm{max}}$. *

Theorem 2.11** (Equivalence of Strong Connectivity and Irreducibility [4]).**

*Consider the (possibly weighted and directed) network that is associated with a nonnegative square matrix $\mathbf{C}$. (That is, $\mathbf{C}$ is the adjacency matrix of the associated network.) The matrix $\bf{C}$ is irreducible if and only if the associated network is strongly connected (i.e., if and only if there exists a path from any origin node to any destination node). *

One typically seeks a centrality matrix that is irreducible (or, equivalently, a matrix for which the associated network that defined by the weighted edges $\{(i,j,C_{ij})\}$ is strongly connected). Ensuring irreducibility is an important consideration when introducing new types of centrality (including ones with both positive and negative edges [28]). For example, the term (1-\alpha)N^{-1}\mathbf{1}\mathbf{1}^{{\color[rgb]{0,0,0}*}} in Definition 2.6 implies that ${\bf C}^{(PR)}$ is positive (i.e., ${C}^{(PR)}_{ij}>0$ for every $i,j\in\mathcal{V}$), which ensures that ${\bf C}^{(PR)}$ is irreducible, regardless of whether the network with adjacency matrix $\mathbf{A}$ is strongly connected [33].

Before continuing, we highlight an eigenvector-based centrality measure that uses both right and left eigenvectors and therefore does not exactly fit Definition 2.2. One defines the so-called dynamical importance of a node in terms of the change of the leading eigenvalue of ${\bf A}$ under removal of that node from the associated network [92] (see also [115]). In practice, as shown in [92], one can approximate dynamical importance to first order (provided one does not lose strong connectivity when removing the node) with an expression that depends on the leading right and left eigenvectors of ${\bf A}$. Other eigenvector-based centralities that involve two or more eigenvectors that one obtains through matrix perturbations have been developed to cater to particular applications, including the spread of infectious diseases [112, 109], percolation [92, 113], and synchronization [114]. Such analysis of perturbations of dynamical systems on networks is also related to notions of eigenvalue and eigenvector elasticities [50, 51, 36].

2.2 Multiplex and Temporal Networks

The different layers of a multilayer network can encode different types of connections and/or interacting systems [55], including interconnected infrastructures [41], categories of social ties [59], networks at different instances of time [118], and many others. By considering the various possibilities for interactions between nodes within and across layers, one can obtain a taxonomy for different types of multilayer networks [55]. We focus on two popular situations: multiplex networks, in which different layers represent different types of interactions; and temporal networks, in which layers represent different time instances or time periods.

We give formal definitions that are salient for multiplex and temporal networks. For both types of multilayer networks, it is convenient to refer to a node $i$ in a layer $t$ as a node-layer pair $(i,t)$.

Definition 2.12** (Uniformly and Diagonally Coupled Multiplex Network).**

*Let $\mathcal{G}(\mathcal{V},\{\mathcal{E}^{(t)}\},\tilde{\mathcal{E}})$ be a multilayer network with nodes $\mathcal{V}=\{1,\dots,N\}$ and $T$ layers, with interactions between node-layer pairs encoded by the sets $\{\mathcal{E}^{(t)}\}$ of weighted edges, where $(i,j,w^{t}_{ij})\in\mathcal{E}^{(t)}$ if and only if there is an edge $(i,j)$ with weight $w^{t}_{ij}$ in layer $t$. The set $\tilde{\mathcal{E}}=\{(s,t,\tilde{w}_{st})\}$ encodes the topology and weights for coupling separate instances of the same node between a layer pair $(s,t)\in\{1,\dots,T\}\times\{1,\dots,T\}$. Equivalently, one can encode a multiplex network as a set $\{\mathbf{A}^{(t)}\}$ of adjacency matrices, such that ${A}_{ij}^{(t)}=w^{t}_{ij}$ if $(i,j,w^{t}_{ij})\in\mathcal{E}^{(t)}$ and ${A}_{ij}^{(t)}=0$ otherwise, along with an interlayer-adjacency matrix $\tilde{\bf{A}}$ with components $\tilde{A}_{st}=\tilde{w}_{st}$ if $(s,t,\tilde{w}_{st})\in\tilde{\mathcal{E}}$ and $\tilde{A}_{ij}^{(t)}=0$ otherwise. *

See Fig. 1(a) for a pedagogical example of a small multiplex network. The multiplex coupling in Definition 2.12 is “diagonal” because we only allow coupling between a node in one layer and that same node in other layers, and it is “uniform” because the coupling between two layers is identical for all nodes in those two layers. A multilayer network with both of these conditions is called “layer-coupled” [55]. Our choice to represent interlayer couplings via $\tilde{\bf{A}}$ is a generalization of the special, but common, case in which the interlayer edge weights are identical for all layer pairs (i.e., $\tilde{\omega}_{st}=\omega$ for all $s$ and $t$). Although there are many other coupling strategies [9, 55], we focus on uniform and diagonal coupling because it is one of the simplest and most popular coupling schemes. The interlayer-adjacency matrix $\tilde{\bf{A}}$ already allows a great deal of flexibility, and (as we will describe in Sec. 3) these restrictions impose matrix symmetries that we can exploit to derive results when the layers are coupled either very weakly or very strongly.

We use a similar multilayer network representation to study temporal networks.

Definition 2.13** (Discrete-Time Temporal Network).**

*A discrete-time temporal network consists of nodes $\mathcal{V}=\{1,\dots,N\}$ and a sequence of network layers. We denote the network either by $\mathcal{G}(\mathcal{V},\{\mathcal{E}^{(t)}\})$ or by a sequence $\{\mathbf{A}^{(t)}\}$ of adjacency matrices such that ${A}_{ij}^{(t)}=w^{t}_{ij}$ if $(i,j,w^{t}_{ij})\in\mathcal{E}^{(t)}$ and ${A}_{ij}^{(t)}=0$ otherwise. *

Definition 2.13 makes no explicit assumptions about how the layers are coupled. That is, a discrete-time temporal network consists of a set of nodes and an ordered set of layers. It is common, however, for temporal networks to also include coupling between layers, such as by representing them (as we do) as a multiplex network with a diagonal interlayer coupling that respects the arrow of time.

2.3 Extensions of Centrality for Multiplex and Temporal Networks

There has been a recent explosion of research on centrality measures for multilayer networks [55, 18]. Much of this work is related to work on generalizing network properties such as node degree [68, 21, 7, 105] and shortest paths, with the latter leading to generalizations of notions like betweenness centrality [67, 99, 100, 14] and closeness centrality [68, 100]. Of particular relevance to the present paper are generalizations of eigenvector-based centralities to multiplex networks. Salient notions that have been generalized include eigenvector centrality [21, 98, 7, 24, 23], hub and authority scores [58, 102, 90, 116], and PageRank [79, 42, 25]. These extensions have employed various strategies; we briefly discuss several of them.

One strategy is to represent a multiplex network as a tensor and use tensor decompositions [58, 122]. Another strategy is to define a system of centrality dependencies in which large-centrality elements (nodes, layers, and so on) connect to other large-centrality elements, and one simultaneously solves for multiple types of centrality [90, 102, 116] as a fixed-point solution of the system of (possibly nonlinear) dependencies. For example, [90] and [116] defined centralities for both nodes and layers such that highly-ranked layers have highly-ranked nodes and highly-ranked nodes are in highly-ranked layers. In a third strategy, which is the one that most closely resembles our present approach, one constructs a supramatrix of size $NT\times NT$ for $N$ nodes and $T$ layers, such that the dominant eigenvector of the supramatrix gives the centralities of the node-layer pairs $\{(i,t)\}$. For example, [98, 94, 23, 24] constructed a supramatrix by using the Khatri–Rao matrix product between a matrix that encodes interlayer connections and a block matrix that has the layers’ adjacency matrices as its block columns. Another approach involves computing one or more centralities independently for each layer and then using a consensus ranking [86]. One can also generalize monolayer centrality measures to construct multilayer “versatility” measures [22].

There have also been many efforts to generalize centralities to temporal networks. There are extensive discussions of such efforts in [110] and [65]. To add to these lists, we briefly highlight several contributions that appeared recently or were not mentioned in [110]. Arrigo and Higham [3] introduced a method to efficiently estimate temporal communicability (a generalization of Katz centrality), which has been applied to a variety of applications, including in neuroscience [69] and the spread of infectious diseases [70, 16]. Huang and Yu [46] extended a measure called dynamic-sensitive centrality to temporal networks. References [87, 47, 30] introduced variants of eigenvector centrality for temporal networks. We highlight [30] in particular, because it explored connections between continuous and discrete-time calculations of temporal centralities. Nathan et al. [76] introduced an efficient algorithm for computing a centrality for streaming graphs. Although methods for streaming and continuous-time networks are important (see, e.g., [1] for a generalization of PageRank to such situations), we restrict our attention to discrete-time temporal networks, especially because we seek to further bridge the literatures on temporal and multiplex networks.

In Sec. 3, we introduce a new construction, which is based on a Kronecker product, for a supracentrality matrix. This construction generalizes our previous work [110], where we introduced a supracentrality framework for temporal networks and assumed adjacent-in-time coupling (i.e., that $\tilde{A}_{tt^{\prime}}=1$ for $|t-t^{\prime}|=1$ and $\tilde{A}_{tt^{\prime}}=0$ otherwise). Our new formulation introduces an interlayer-adjacency matrix $\tilde{\bf{A}}$, allowing our framework to flexibly cater to either multiplex or temporal networks. Note that there are also multilayer representations of temporal networks that are not multiplex. One such approach is to connect node-layer pair $(i,t)$ to $\{(j,t+1)\}$ for $j\in\{j:A^{(t)}_{ij}\not=0\}$, yielding a supra-adjacency matrix with identity matrices on the block diagonal and the layers’ adjacency matrices on the off-diagonal blocks that lie directly above the diagonal blocks. These interlayer edges are nondiagonal, because they connect nodes in one layer to different nodes in a neighboring layer. This formulation, which is connected mathematically [29] to matrix-iteration-based centrality measures for temporal networks [37, 38, 39], has been used to study time-dependent functional brain networks [117]. Additionally, multilayer networks with non-diagonal edges were used in [118] to study disease spreading on temporal networks. One can also to choose to study multiplex and temporal networks without interlayer coupling by independently considering each layer in isolation.

3 Supracentrality Framework for Multiplex and Temporal Networks

We now present a supracentrality framework that provides a common mathematical foundation for eigenvector-based centralities for layer-coupled multiplex and temporal networks. In Sec. 3.1, we define supracentrality matrices. In Sec. 3.2, we define joint, marginal, and conditional centralities; we prove their uniqueness and positivity under certain conditions. In Sec. 3.3, we give a pedagogical example to illustrate these types of centralities.

We summarize our key mathematical notation in Table 1. We use the subscripts {{\color[rgb]{0,0,0}i,j}}\in\mathcal{V} to enumerate nodes, the subscripts {{\color[rgb]{0,0,0}s,t}}\in\mathcal{\{}1,\dots,T\} to enumerate layers, and the subscripts {{\color[rgb]{0,0,0}p,q}}\in\mathcal{\{}1,\dots,NT\} to enumerate node-layer pairs.

3.1 Supracentrality Matrices

We first define a supracentrality matrix, in a way that generalizes the definition in [110], for networks that are either multiplex or temporal.

Definition 3.1** (Supracentrality Matrix).**

Let $\{{\bf C}^{(t)}\}$ be a set of $T$ centrality matrices for a multilayer network whose layers have a common set $\mathcal{V}=\{1,\dots,N\}$ of nodes, and suppose that ${C}^{(t)}_{ij}\geq 0$. Let $\tilde{\bm{A}}$, with $\tilde{A}_{ij}\geq 0$, be a $T\times T$ interlayer-adjacency matrix that encodes the interlayer couplings. We define a family of supracentrality matrices $\mathbb{C}(\omega)$, which are parameterized by the interlayer-coupling strength $\omega\geq 0$, of the form

[TABLE]

*where $\hat{\mathbb{C}}=\text{\rm diag}[\mathbf{C}^{(1)},\dots,\mathbf{C}^{(T)}]$ and $\hat{\mathbb{A}}=\tilde{\bm{A}}\otimes\bf{I}$ denotes the Kronecker product of $\tilde{\bm{A}}$ and $\bf{I}$. *

Remark 3.2**.**

*For layer $t$, the matrix ${\bf C}^{(t)}$ can be any matrix whose dominant eigenvector is of interest. We focus on centrality matrices, such as those that are associated with eigenvector centrality (see Definition 2.3), hub and authority scores (see Definition 2.4), and PageRank (see Definition 2.6). Additionally, one can scale each ${\bf C}^{(t)}$ by a layer-specific weight. (Such a scaling has the potential to benefit both multilayer community detection [81] and layer-averaged clique detection [77].) One can easily incorporate such weighting into Eq. (11) by redefining the centrality matrices $\{{\bf C}^{(t)}\}$ to include the weights. *

The supracentrality matrix $\mathbb{C}(\omega)$ of size $NT\times NT$ encodes the effects of two distinct types of connections: the layer-specific centrality entries $\{{C}^{(t)}_{ij}\}$ in the diagonal blocks relate centralities between nodes within layer $t$, whereas entries in the off-diagonal blocks encode coupling between layers. The supramatrix $\hat{\mathbb{A}}=\tilde{\bm{A}}\otimes\bf{I}$ implements uniform and diagonal coupling: the matrix $\bf{I}$ encodes diagonal coupling, and any two layers $t$ and $t^{\prime}$ are uniformly coupled, because all interlayer edges between them have the identical weight $w\tilde{A}_{tt^{\prime}}$. The choice of undirected, adjacent-in-time interlayer coupling (i.e., $\tilde{A}_{tt^{\prime}}=1$ if $|t-t^{\prime}|=1$ and $\tilde{A}_{tt^{\prime}}=0$ otherwise) recovers the supracentrality matrix that we studied in [110]. In the present paper, we generalize that notion of a supracentrality matrix by using an interlayer-adjacency matrix $\tilde{\bm{A}}$, which allows us to implement a wide variety of interlayer coupling topologies. In the context of multiplex networks, we hypothesize that different choices for $\tilde{\bm{A}}$ will have different benefits. In the context of temporal networks, we will study (see Sec. 5.2) the effects of letting $\tilde{\bm{A}}$ encode a directed, time-respecting chain with “layer teleportation” (see Eq. (29)). This yields supracentrality results that we will contrast with those in [110].

3.2 Joint, Marginal, and Conditional Centralities

We now study supracentralities in the form of joint, marginal, and conditional centralities.

The defining feature of eigenvector-based centrality measures is that one computes and studies a dominant eigenvector of a centrality matrix. We study the right dominant-eigenvalue equation

[TABLE]

where we interpret entries in the right dominant eigenvector $\mathbbm{v}(\omega)$ as centrality measures for the node-layer pairs $\{(i,t)\}$. The vector $\mathbbm{v}(\omega)$ has a block form: its first $N$ entries encode the joint centralities for layer $t=1$, its next $N$ entries encode the joint centralities for layer $t=2$, and so on. Therefore, as we now describe, it can be useful to reshape the block vector $\mathbbm{v}(\omega)$ into a matrix.

Following [110], we use the concepts of joint, marginal, and conditional centralities to develop our understanding of the importances of nodes and layers from the values of $\mathbbm{v}(\omega)$.

Definition 3.3** (Joint Centralities of Node-Layer Pairs [110]).**

Let $\mathbb{C}(\omega)$ be a supracentrality matrix given by Definition 3.1, and let $\mathbbm{v}(\omega)$ be its right dominant eigenvector. We encode the joint centrality of node $i$ in layer $t$ via the $N\times T$ matrix ${\bf W}(\omega)$ with entries

[TABLE]

Remark 3.4**.**

*We refer to $W_{it}(\omega)$ as a “joint centrality” because it reflects the importance both of node $i$ and of layer $t$. *

Definition 3.5** (Marginal Centralities of Nodes and Layers [110]).**

Let ${\bf W}(\omega)$ encode the joint centralities given by Definition 3.3. We define the marginal layer centrality (MLC) ${x}_{t}(\omega)$ and marginal node centrality (MNC) $\hat{x}_{i}(\omega)$ by

[TABLE]

Definition 3.6** (Conditional Centralities of Nodes and Layers [110]).**

Let the set $\{W_{it}(\omega)\}$ be the joint centralities given by Definition 3.3, and let $\{{x}_{t}(\omega)\}$ and $\{\hat{x}_{i}(\omega)\}$, respectively, be the marginal layer and node centralities from Definition 3.5. We define the conditional centralities of nodes and layers by

[TABLE]

*where $Z_{it}(\omega)$ gives the centrality of node $i$ conditioned on layer $t$ and $\hat{Z}_{it}(\omega)$ gives the centrality of layer $t$ conditioned on node $i$. *

The quantity $Z_{it}(\omega)$ indicates the importance of node $i$ relative to other nodes in layer $t$. By contrast, the joint node-layer centrality $W_{it}(\omega)$ measures the importance of node-layer pair $(i,t)$ relative to all node-layer pairs.

We now present a new theorem that ensures the uniqueness and positivity of the above supracentralities.

Theorem 3.7** (Uniqueness and Positivity of Supracentralities).**

Let $\mathbb{C}(\omega)$ be a supracentrality matrix given by Eq. (11). Additionally, suppose that $\tilde{\bm{A}}$ is an adjacency matrix for a strongly connected graph and that the sum $\sum_{t}\mathbf{C}^{(t)}$ is an irreducible nonnegative matrix. It then follows that $\mathbb{C}(\omega)$ is irreducible, nonnegative, and has a simple largest positive eigenvalue $\lambda_{\rm{max}}(\omega)$ with corresponding left eigenvector $\mathbbm{u}(\omega)$ and right eigenvector $\mathbbm{v}(\omega)$, which are unique and consist of positive entries. Moreover, the centralities $\{W_{it}(\omega)\}$, $\{x_{t}(\omega)\}$, $\{\hat{x}_{i}(\omega)\}$, $\{Z_{it}(\omega)\}$, and $\{\hat{Z}_{it}(\omega)\}$ are positive and well-defined.

Proof 3.8**.**

*See Appendix A. *

Remark 3.9**.**

*If we also assume that $\mathbb{C}(\omega)$ is aperiodic, then it is also true that $\lambda_{\rm{max}}(\omega)$ is larger in magnitude than the other eigenvalues. *

3.3 Pedagogical Example that Illustrates Different Coupling Regimes

In Fig. 2, we illustrate the concepts of joint and marginal centralities for the multiplex network in Fig. 1(a). This network has $N=4$ nodes and $T=6$ layers, and we study the interlayer-adjacency matrix $\tilde{\bm{A}}$ that we showed in the inset of Fig. 1(a). We set $\tilde{A}_{tt^{\prime}}=1$ for all depicted interlayer couplings, except for the coupling of layers $3$ and $4$, for which we set $\tilde{A}_{34}=\tilde{A}_{43}=0.01$. With these interlayer edge weights, the interlayer-coupling network that is associated with $\tilde{\bm{A}}$ has two natural communities of densely-connected nodes. For this experiment (and our other experiments), we typically find that conditional node centralities provide the most useful insights.

In Fig. 3(a), we plot the conditional centralities $\{Z_{it}(\omega)\}$ of node-layers for three different choices of the interlayer coupling-strength $\omega$. These choices represent three centrality regimes (which we illustrate in Fig. 3(b)) that we observe by exploring centralities across a range of $\omega$ values. In the top two panels of Fig. 3(b), we plot the MNC and MLC values versus $\omega$. In the bottom panel, we quantify the “sensitivity” of the joint and conditional centralities to perturbations of $\omega$. Specifically, we consider $\omega$ in the interval $[10^{-2},10^{4}]$ discretized by $\omega_{s}=10^{-2+0.2s}$ for $s\in\{0,\dots,30\}$. We plot the stepwise magnitudes of the changes $\|{\bf W}(\omega_{s})-{\bf W}(\omega_{s-1})\|_{F}$ of the joint centralities and the changes $\|{\bf Z}(\omega_{s})-{\bf Z}(\omega_{s-1})\|_{F}$ of the conditional centralities, where $\|\cdot\|_{F}$ denotes the Frobenius norm. We identify three regimes for which the conditional centralities are robust. (See the shaded regions in the bottom panel of Fig. 3(b).) The peaks in Fig. 3(b) indicate values of $\omega$ where conditional centralities are most sensitive to perturbations of $\omega$; other choices for $\omega$ are more robust to a perturbation of $\omega$. Interestingly, the peaks and troughs for the curves for $\|{\bf W}(\omega_{s})-{\bf W}(\omega_{s-1})\|_{F}$ and $\|{\bf Z}(\omega_{s})-{\bf Z}(\omega_{s-1})\|_{F}$ do not coincide. We focus on robust values of ${\bf Z}(\omega)$, because we generally find that conditional centralities provide the most interpretable and insightful results among our supracentralities. See Sec. 3.2 and Sec. 5 of [110].

We summarize these three regimes as follows.

• •

Weak-Coupling Regime. The top panel of Fig. 3(a), which shows results for interlayer-coupling strength $\omega=10^{-2}$, depicts a regime in which the conditional node centralities resemble the centralities of uncoupled layers as $\omega\rightarrow 0^{+}$. In this example, we observe that the conditional centralities correlate strongly with intralayer degrees ($d_{i}^{(t)}=\sum_{j}A^{(t)}_{ij}$). Specifically, as we indicate with the horizontal dashed lines, the conditional centralities approximately equal one of three values: $Z_{it}\approx 0.3656$ for $d_{i}^{(t)}=3$, $Z_{it}\approx 0.3096$ for $d_{i}^{(t)}=2$, and $Z_{it}\approx 0.2021$ for $d_{i}^{(t)}=1$. We analyze the $\omega\rightarrow 0^{+}$ limit in Sec. 4.1.

• •

Strong-Coupling Regime. The bottom panel of Fig. 3(a), which shows results for $\omega=10^{4}$, depicts a regime in which the conditional centralities approach the centralities of a layer-aggregated centrality matrix as $\omega\rightarrow\infty$. In this regime, the conditional centralities $Z_{it}$ of each node $i$ limit to a value $\alpha_{i}$ that is constant across the layers; that is, $Z_{it}\to\alpha_{i}$ for all $t$. We analyze the $\omega\rightarrow\infty$ limit in Sec. 4.2.

• •

Intermediate-Coupling Regime. The middle panel of Fig. 3(a) depicts an intermediate regime. We show conditional centralities at $\omega=10$, which is a value of $\omega$ that lies between the locations of the two peaks of the bimodal curve $\|{\bf Z}(\omega_{s})-{\bf Z}(\omega_{s-1})\|_{F}$ that we show in Fig. 3(b). For $\omega=10$, we observe two subsets of the $Z_{it}$ conditional centrality values: those for layers $t\in\{1,2,3\}$ are very similar to each other, and those for layers $t\in\{4,5,6\}$ are also similar to each other. This pattern arises directly from the layer-coupling scheme in Fig. 1(a), where these two sets of layers correspond to two communities in the interlayer-adjacency matrix. (Recall that the coupling between layers $3$ and $4$ is $100$ times weaker than the other couplings.) In Sec. SM1 of the Supplementary Materials, we show that the curve $\|{\bf Z}(\omega_{s})-{\bf Z}(\omega_{s-1})\|_{F}$ becomes unimodal as we increase the coupling between layers $3$ and $4$.

We classify the strong-coupling and weak-coupling regimes by considering whether or not the observed supracentralities are strongly correlated with those of either asymptotic limit (i.e., either $\omega\to 0^{+}$ or $\omega\to\infty$). The intermediate regime arises from an interplay between the topologies and edge weights of the layers and the interlayer couplings, so these multilayer centralities provide insights that cannot be observed by studying the network layers in isolation or in aggregate. This example also illustrates that it is important to explore various coupling strengths $\omega$ and various interlayer-adjacency matrices $\tilde{\bm{A}}$ to identify supracentralities that are appropriate for a given application. See [106] for our Matlab code that computes supracentralities and reproduces Fig. 3 and our other experiments in this paper.

4 Limiting Behavior for Weak and Strong Coupling

We construct singular perturbation expansions to analyze the limiting behaviors of Eq. (12) when the interlayer-coupling strength $\omega$ is very small (i.e., layer decoupling) or very large (i.e., layer aggregation). These results provide insights into our supracentrality framework and can aid in the selection of appropriate parameter values.

4.1 Layer Decoupling in the Weak-Coupling Limit

We first study supracentralities in the $\omega\to 0^{+}$ limit. We did not study this limit in our previous work [110] on temporal centralities. We analyze how the eigenvalues and eigenvectors of $\mathbb{C}(\omega)$ for small $\omega$ are determined by the eigenvalues and eigenvectors of the individual layers’ centrality matrices ${\bf C}^{(t)}$.

Lemma 4.1** (Layer Decoupling at $\omega=0$).**

Let $\mu_{i}^{(t)}$, ${\bf v}^{(i,t)}$ and ${\bf u}^{(i,t)}$ denote the eigenvalues and corresponding right and left eigenvectors of the $N\times N$ centrality matrices ${\bf C}^{(t)}$ for $t\in\{1,\dots,T\}$. It follows that each $\mu_{i}^{(t)}$ is an eigenvalue of the associated supracentrality matrix $\mathbb{C}(0)=\mathrm{diag}[{\bf C}^{(1)}$,…,${\bf C}^{(T)}]$ with corresponding left eigenvector $\mathbbm{u}^{(i,t)}={\bf e}^{(t)}\otimes{\bf u}^{(i,t)}$ and right eigenvector $\mathbbm{v}^{(i,t)}={\bf e}^{(t)}\otimes{\bf v}^{(i,t)}$, where ${\bf e}^{(t)}$ denotes a length-$T$ unit vector that consists of zeros in all entries except for the $t$-th entry (which is a $1$).

Proof 4.2**.**

*See Appendix B. *

Remark 4.3**.**

*Each of the vectors $\mathbbm{u}^{(i,t)}$ and $\mathbbm{v}^{(i,t)}$ consists of $T$ blocks, and each block is a length-$N$ vector that consists of zeros except for the $t$-th block, which is the associated left or right eigenvector of ${\bf C}^{(t)}$. *

Remark 4.4**.**

If there exist multiple eigenvalues $\mu_{i}^{(t)}$ of ${\bf C}^{(t)}$ that are equal (that is, $\mu_{i}^{(t)}=\lambda$ for any $(i,t)$ in some set $\mathcal{P}$), then any linear combination of their associated eigenvectors is also an eigenvector:

[TABLE]

*One can show a similar result for the left eigenvectors. Consequently, the eigenspace that is associated with each eigenvalue is an invariant subspace. *

We now turn our attention to the eigenspace of the dominant eigenvalue (i.e., the dominant eigenspace). Let $\lambda_{\rm{max}}(0)$ denote the largest positive (i.e., dominant) eigenvalue of $\mathbb{C}(0)$, and let $\mathcal{P}=\{t:\mu_{1}^{(t)}=\lambda_{\rm{max}}(0)\}$ denote the set of indices for layers whose largest positive eigenvalue $\mu_{1}^{(t)}$ is equal to $\lambda_{\rm{max}}(0)$. We assume that each layer’s dominant eigenvalue $\mu_{1}^{(t)}$ is simple (i.e, $\mu_{j}^{(t)}<\mu_{1}^{(t)}$ for any $j\geq 2$) and that its right and left eigenvectors are unique. This implies that $i=j=1$ for our notation in Eq. (4.4), so we can consider a set $\mathcal{T}$ of layers rather than a set $\mathcal{P}$ of node-layer pairs.

We now present a key result for the dominant eigenvectors of $\lim_{\omega\to 0^{+}}\mathbb{C}(\omega)$.

Theorem 4.5** (Weak-Coupling Limit of Dominant Eigenvectors).**

Let $\mathbbm{v}^{(1)}(\omega)$ and $\mathbbm{u}^{(1)}(\omega)$, respectively, be the right and left dominant eigenvectors of a supracentrality matrix $\mathbb{C}(\omega)$ under the assumptions of Thm. 3.7. It then follows that the $\omega\to 0^{+}$ limits of $\mathbbm{u}^{(1)}(\omega)$ and $\mathbbm{v}^{(1)}(\omega)$ are

[TABLE]

where the vectors \bm{\alpha}=[\alpha_{1},\dots,\alpha_{T}]^{{\color[rgb]{0,0,0}*}} and \bm{\beta}=[{\beta}_{1},\dots,\beta_{T}]^{{\color[rgb]{0,0,0}*}}, which have nonnegative entries that satisfy $\sum_{t}\alpha_{t}^{2}=\sum_{t}\beta_{t}^{2}=1$, are positive solutions to the dominant-eigenvalue equations

[TABLE]

where $\lambda_{1}$ is an eigenvalue that needs to be determined, the entries of $\bm{X}$ are

[TABLE]

and \chi(t)=\sum_{t^{\prime}\in{{\color[rgb]{0,0,0}\mathcal{T}}}}\delta_{tt^{\prime}} is an indicator function (with $\chi(t)=1$ if t\in{{\color[rgb]{0,0,0}\mathcal{T}}} and $\chi(t)=0$ otherwise). (Recall that $*$ is our notation for transpose operator.)

Proof 4.6**.**

*See Appendix C. *

Remark 4.7**.**

*We obtained Eq. (19) using singular perturbation theory in the limit $\omega\to 0^{+}$. One can understand why the limit is singular by considering the dimension $f(\omega)=\mathrm{dim}(\mathrm{null}(\mathbb{C}(\omega)-\lambda_{\rm{max}}(\omega)\mathbb{I}))$ of the eigenspace that is associated with the dominant eigenvalue $\lambda_{\rm{max}}(\omega)$ of $\mathbb{C}(\omega)$. For any $\omega>0$, Thm. 3.7 (i.e., our Perron–Frobenius theorem) guarantees that $f(\omega)=1$. However, Eq. (4.4) implies that when $\omega=0$, there are |{{\color[rgb]{0,0,0}\mathcal{T}}}| eigenvectors that are associated with eigenvalue $\lambda_{\rm{max}}(0)$, which in turn implies that f(0)=|{{\color[rgb]{0,0,0}\mathcal{T}}}|. If |{{\color[rgb]{0,0,0}\mathcal{T}}}|>1, then $\lim_{\omega\to 0^{+}}f(\omega)\not=f(0)$ and $\omega=0$ is a singular point of the dominant eigenspace. *

We now present three corollaries that consider Thm. 4.5 under various restrictions on the centrality matrices. We first consider the limiting behavior when the layers’ centrality matrices all have the same spectral radius, as is the case for PageRank matrices (because a PageRank matrix is a transition matrix of a Markov chain) or if one rescales the centrality matrices to have the same spectral radius.

Corollary 4.8** (Weak-Coupling Limit for

Centrality Matrices with the Same Spectral Radius).**

Under the assumptions of Thm. 4.5 and also assuming that all centrality matrices have the same spectral radius (i.e., $\lambda_{\rm{max}}=\mu_{1}^{(t)}$ for all $t$), it follows that {{\color[rgb]{0,0,0}\mathcal{T}}}=\{1,\dots,T\} and $\chi(t)=1$. Additionally, Eq. (19) takes the form

[TABLE]

We next consider when the layers’ centrality matrices are symmetric, which is the case for hub/authority matrices and for symmetric adjacency matrices.

Corollary 4.9** (Weak-Coupling Limit for Symmetric Centrality Matrices).**

Under the assumptions of Thm. 4.5 and also assuming that all centrality matrices are symmetric, ${\bf u}^{(1,t)}={\bf v}^{(1,t^{\prime})}$ and Eq. (19) takes the form

[TABLE]

When the centrality matrix of a single layer has the largest spectral radius, which one often expects to occur for adjacency matrices and hub/authority matrices (unless the network layers have symmetries that yield repeated spectral radii across layers), the limiting behavior of the eigenvector is that it localizes onto a single dominating layer.

Corollary 4.10** (Weak-Coupling-Induced Eigenvector Localization onto a Dominating Layer).**

Under the assumptions of Thm. 4.5 and also assuming that one layer has a spectral radius that is larger than all others (i.e., $\lambda_{\rm{max}}=\mu_{1}^{(t)}$ for a single layer t={{\color[rgb]{0,0,0}\tau}}), then as $\omega\to 0^{+}$, it follows that

[TABLE]

Understanding whether or not the dominant eigenvector localizes onto a single layer, onto several layers (i.e., as given by the function $\chi(t)$), or does not localize has significant practical consequences. In some situations, it can be appropriate to allow eigenvector localization [83], whereas it can be beneficial to avoid localization in others [71]. Lemma 4.1 and Corollaries 4.9–4.10 characterize localization in the weak-coupling limit and are useful for practitioners to make informed choices about which centrality matrices to use.

4.2 Layer Aggregation in the Strong-Coupling Limit

We study Eq. (12) in the limit as $\omega\to\infty$ (or, equivalently, as $\epsilon:=\omega^{-1}\to 0^{+}$). The results of this subsection generalize those of [110], where we assumed that $\tilde{\bm{A}}$ encodes adjacent-in-time coupling. In the present discussion, by contrast, we allow $\tilde{\bm{A}}$ to be from a much more general class of matrices, including asymmetric matrices that encode directed interlayer couplings.

Consider the scaled supracentrality matrix

[TABLE]

which has eigenvectors $\tilde{\mathbbm{u}}(\epsilon)$ and $\tilde{\mathbbm{v}}(\epsilon)$ that are identical to those of $\mathbb{C}(\omega)$ (specifically, $\tilde{\mathbbm{u}}(\epsilon)={\mathbbm{u}}(\epsilon^{-1})$ and $\tilde{\mathbbm{v}}(\epsilon)={\mathbbm{v}}(\epsilon^{-1})$). Its eigenvalues $\{\tilde{\lambda}_{i}\}$ are scaled by $\epsilon$; specifically, $\tilde{\lambda}_{i}(\epsilon)=\epsilon\lambda_{i}(\epsilon^{-1})$.

To facilitate our presentation, we define a permutation operator for $NT\times NT$ matrices.

Definition 4.11** (Node-Layer-Reordering Stride Permutation).**

The matrix  $\mathbb{P}$ is a $T$-stride permutation matrix of size $NT\times NT$ if it has entries that take the form [34]

[TABLE]

*Therefore, \left(\tilde{\bm{A}}\otimes\mathbf{I}\right)=\mathbb{P}\left(\mathbf{I}\otimes\tilde{\bm{A}}\right)\mathbb{P}^{{\color[rgb]{0,0,0}*}}. *

Remark 4.12**.**

*The stride-permutation matrix is unitary, and it simply changes the ordering of node-layer pairs. Before the stride permutation (which is a type of graph isomorphism [56]), a supracentrality matrix has entries that are ordered first by node $i$ and then by layer $t$ (i.e., $(i,t)=(1,1),(2,1),(3,1),\dots$). After the stride permutation, the entries are ordered first by layer $t$ and then by node $i$ (i.e., $(i,t)=(1,1),(1,2),(1,3),\dots$). *

We now present our main findings for the strong-coupling regime.

Lemma 4.13** (Singularity at Infinite Coupling).**

Let $\tilde{\mu}_{t}$ denote the eigenvalues of $\tilde{\bm{A}}$, and let $\tilde{\bm{v}}^{(t)}$ and $\tilde{\bm{u}}^{(t)}$, respectively, be the corresponding right and left eigenvectors. We assume that the eigenvalues are simple, and we order them such that $\tilde{\mu}_{1}$ is the largest eigenvalue. We also let $\mathbb{P}$ denote the stride permutation matrix from Eq. (24).

For $\epsilon=0$, each $\tilde{\mu}_{t}$ is an eigenvalue of $\mathbb{C}(\epsilon)$ given by Eq. (23), and the associated $N$-dimensional right and left eigenspaces are spanned by the eigenvectors $\mathbb{P}\tilde{\mathbbm{v}}^{(t)}$ and $\mathbb{P}\tilde{\mathbbm{u}}^{(t)}$, respectively, with the general form

[TABLE]

where the constants $\tilde{\alpha}_{tj}$ and $\tilde{\beta}_{tj}$ must satisfy $\sum_{j}\tilde{\alpha}_{tj}^{2}=\sum_{j}\tilde{\beta}_{tj}^{2}=1$ to ensure that $\|\tilde{\mathbbm{u}}^{(t)}\|_{2}=\|\tilde{\mathbbm{v}}^{(t)}\|_{2}=1$. The associated length-$NT$ vectors are $\tilde{\mathbbm{v}}^{(t,j)}=\tilde{{\bf e}}^{(j)}\otimes{\bf v}^{(t)}$ and $\tilde{\mathbbm{u}}^{(t,j)}=\tilde{{\bf e}}^{(j)}\otimes{\bf u}^{(t)}$, where $\tilde{{\bf e}}^{(j)}$ is a length-$N$ unit vector that consists of zeros in all entries except for entry $j$, which is $1$. Therefore, $\tilde{\mathbbm{u}}^{(t,j)}$ (respectively, $\tilde{\mathbbm{v}}^{(t,j)}$) consists of zeros, except in the $j$-th block of size $T$, which consists of a left (respectively, right) eigenvector of $\tilde{\bm{A}}$.

Proof 4.14**.**

*See Appendix D. *

Remark 4.15**.**

*It is straightforward to also obtain the general form of eigenvectors for eigenvalues $\{\mu_{t}\}$ whose multiplicity is larger than $1$. For example, if the eigenvalue $\tilde{\mu}_{t}$ of $\tilde{\bm{A}}$ has multiplicity $q$, then $\tilde{\lambda}_{i}(0)=\tilde{\mu}_{t}$ has multiplicity $qN$ for the matrix $\mathbb{C}(0)$. However, the notation becomes slightly more cumbersome, and we will not study such cases in this paper. *

Theorem 4.16** (Strong-Coupling Limit of Dominant Eigenvectors).**

Let $\tilde{\mu}_{1}$ denote the dominant eigenvalue (which we assume to be simple) of the interlayer-adjacency matrix $\tilde{\bm{A}}$, and let $\tilde{\bm{v}}^{(1)}$ and $\tilde{\bm{u}}^{(1)}$ be its associated right and left eigenvectors. We assume that the constraints of Thm. 3.7 are satisfied, such that the supracentrality matrix $\mathbb{C}(\epsilon)$ given by Eq. (11) is nonnegative, irreducible, and aperiodic. It then follows that the largest positive eigenvalue $\tilde{\lambda}_{\rm{max}}(\epsilon)$ and its associated eigenvectors, $\mathbbm{u}^{(1)}(\epsilon)$ and $\mathbbm{v}^{(1)}(\epsilon)$, of $\mathbb{C}(\epsilon)$ converge as $\epsilon\to 0^{+}$ as follows:

[TABLE]

where $\mathbb{P}$ is the stride permutation from Eq. (24), the vectors $\tilde{\mathbbm{u}}^{(1,j)}$ and $\tilde{\mathbbm{v}}^{(1,j)}$ are defined in Thm. 4.13, and the constants $\tilde{\beta}_{i}$ and $\tilde{\alpha}_{i}$ solve the dominant-eigenvalue equations

[TABLE]

where

[TABLE]

Proof 4.17**.**

*See Appendix E. *

Equation (28) indicates that the strong-coupling limit effectively aggregates the centrality matrices $\{{\bf C}^{(t)}\}$ across time via a weighted average, with weights that depend on the right and left dominant eigenvectors of the interlayer-adjacency matrix $\tilde{\bm{A}}$. This result generalizes Eq. (4.13) of Taylor et al. [110], who assumed that $\tilde{\bm{A}}$ is symmetric (so that $\tilde{u}_{t}^{(1)}=\tilde{v}_{t}^{(1)}$). We recover the result in [110] with the following corollary.

Corollary 4.18** (Strong-Coupling Limit of Eigenvector-Based Centralities for Multilayer Networks with Adjacent-in-Time, Uniform, and Diagonal Coupling [110]).**

For undirected, adjacent-in-time interlayer coupling (i.e., $\tilde{A}_{tt^{\prime}}=1$ for $|t-t^{\prime}|=1$ and $\tilde{A}_{tt^{\prime}}=0$ otherwise), the $\epsilon\to 0^{+}$ limit of the largest eigenvalue is $\tilde{\lambda}_{1}(\epsilon)\to 2\cos\left({\frac{\pi}{T+1}}\right)$. When $\epsilon\to 0^{+}$, the right and left dominant eigenvectors satisfy Eqs. (26)–(27), with

[TABLE]

Corollary 4.19** (Strong-Coupling Limit of Eigenvector-Based Centralities for Multilayer Networks with All-to-All, Uniform, and Diagonal Coupling).**

For all-to-all coupling (with self-edges), \tilde{\bm{A}}=\bm{1}\bm{1}^{{\color[rgb]{0,0,0}*}}; the $\epsilon\to 0^{+}$ limit of the largest eigenvalue is $\lambda_{\rm{max}}(\epsilon)\to\tilde{\mu}_{1}=N$; and, when $\epsilon\to 0^{+}$, the right and left dominant eigenvectors satisfy Eqs. (26)–(27), with $\tilde{\mathbf{X}}=T^{-1}\sum_{t}\mathbf{C}^{(t)}$.

Proof 4.20**.**

*In this case, the largest eigenvalue of $\tilde{\bm{A}}$ is $\mu_{1}=N$, and the right and left dominant eigenvectors have the same value in each component, with entries $u_{t}^{(1)}=v_{t}^{(1)}=T^{-1/2}$. *

Corollary 4.21** (Strong-Coupling Limit of of Eigenvector-Based Centralities for Multilayer Networks with Rank-$1$, Uniform, and Diagonal Coupling).**

For a rank-$1$ interlayer-adjacency matrix, \tilde{\bm{A}}=\tilde{\bm{w}}\tilde{\bm{w}}^{{\color[rgb]{0,0,0}*}}, the $\epsilon\to 0^{+}$ limit of the largest eigenvalue is $\lambda_{\rm{max}}(\epsilon)\to 1$; and, when $\epsilon\to 0^{+}$, the right and left dominant eigenvectors satisfy Eqs. (26)–(27), with $\tilde{\mathbf{X}}=\sum_{t}\mathbf{C}^{(t)}\tilde{w}_{t}^{2}$.

Proof 4.22**.**

*This follows from the fact that the largest eigenvalue of $\tilde{\bm{A}}$ is $\tilde{\mu}_{1}=1$, with associated eigenvectors $\tilde{u}_{t}^{(1)}=\tilde{v}_{t}^{(1)}=\tilde{w}_{t}$. *

5 Case Studies with Empirical Data

We now apply our supracentrality framework to study networks that we construct from two data sets. In Sec. 5.1, we examine a multiplex network that encodes flight patterns between European airports, with layers representing different airline companies [13]. In Sec. 5.2, we examine a temporal network that encodes the exchange of mathematicians between mathematical-sciences doctoral programs in the United States [110]. These experiments explore different interlayer-coupling strengths and topologies, and they also illustrate strategies for how to flexibly apply our supracentrality framework to diverse applications.

5.1 Rankings of European Airports in a Multiplex Airline Network [13]

We first apply our supracentrality framework to study the importances of European airports using an empirical multiplex network of the flight patterns for $T=37$ airline companies [13]. See Fig. 4 for a map of the airports and a visualization of the network layers that have the most edges. The network’s nodes represent European airports; we consider only the $N=417$ nodes in the largest connected component of the network that is associated with the sum of the layers’ adjacency matrices. We couple the layers using all-to-all coupling (so \tilde{\bm{A}}=\bm{1}\bm{1}^{{\color[rgb]{0,0,0}*}}), and we calculate eigenvector supracentralities with the centrality matrices from Definition 2.3. See Table 2 for a list of the airline companies that are associated with the network layers. For each layer, we indicate the number $M_{t}=\frac{1}{2}\sum_{ij}A^{(t)}_{ij}$ of intralayer edges and the spectral radius $\lambda_{1}^{(t)}$ of the layer’s associated adjacency matrix ${\bf A}^{(t)}$.

For each airport, we compute the joint, marginal, and conditional centralities for a range of interlayer-coupling strengths $\omega$. In Fig. 5, we plot the MLCs and MNCs (see Definition 3.5). In Fig. 5(a), we see that for large $\omega$, all layers have similar importances; by contrast, for small $\omega$, one layer is much more important because the eigenvector-localization phenomenon that we described in Cor. 4.10. Specifically, the layer that represents Ryanair dominates for small $\omega$, as its adjacency matrix has the largest spectral radius. (It also has the most edges.) A previous investigation of multilayer centralities in this data set [90] also identified Ryanair as the most important airline.

The insets in Fig. 5 give the calculated conditional node centralities for $\omega=10^{-2}$ (in panel (a)) and $\omega=10^{2}$ (in panel (b)) versus the asymptotic values in the associated limits for $\omega\to 0^{+}$ (see Thm. 4.5) and $\omega\to\infty$ (see Thm. 4.16), demonstrating that they are in excellent agreement. One can also observe in Fig. 5(b) that there is not a simple transition between these two limits. Specifically, the thick black curves highlight a few airports whose MNCs have a peak for intermediate values of $\omega$. These airports are more important if one considers the airline network as a multiplex network than if one considers the layers in isolation or in aggregate.

In Table 3, we list the airports with the largest MLCs for eigenvector supracentrality for small, intermediate, and large values of $\omega$. As expected, for large and small $\omega$, the top airports correspond to the top airports (i.e., those with the largest eigenvector centralities) that are associated with the aggregation of layers and the Ryanair network layer, respectively. The top-ranked airports for $\omega=1$ have a large overlap with those for $\omega=100$. The top airport, Adolfo Suárez Madrid–Barajas Airport (LEMD), is particularly interesting. LEMD is the only airport that ranks in the top 10 for both the Ryanair layer and the layer-aggregated network; this contributes to its having the top rank for this intermediate value of $\omega$. We highlight similar ranking boosts for other airports with solid black curves in Fig. 5(b). We also note that LEMD was identified in previous investigations [48, 116] as one of the most important airports in this data set.

In Fig. 6, we illustrate that the eigenvector supracentralities correlate strongly with node degrees. In Fig. 6(a), we show for $\omega=100$ that the airports’ conditional centralities, averaged across layers, are correlated strongly with their total (i.e., layer-aggregated) degrees $\overline{d}_{i}=\sum_{t,j}A_{ij}^{(t)}$ (see the blue ‘$\times$’ marks). We expect this strong correlation for eigenvector supracentrality, as node degree is a first-order approximation of eigenvector centrality in monolayer networks [112]. We also plot the mean conditional centralities versus the number of length-2 paths that emanate from each node (see the red circles). As expected, this correlation is even stronger, as the number $\sum_{t,j}[A^{(t)}]_{ij}^{2}$ of length-2 paths is a second-order approximation to eigenvector centrality [112].444As described in [112], one can interpret the number of length-$k$ paths as an order-$k$ approximation to the dominant eigenvector of an adjacency matrix when the largest positive (i.e., dominant) eigenvalue has a magnitude that is strictly larger than those of all of the other eigenvalues. Recall that Thm. 2.10 guarantees that this true whenever ${\bf A}$ is nonnegative, irreducible, and aperiodic. (We also note that the aperiodic assumption is invalid for certain classes of networks, such as bipartite networks.) Consider the power-method iteration for numerically computing the dominant eigenvector for any adjacency matrix ${\bf A}$ under these assumptions. If one initializes the power method with ${\bf 1}$ (a vector of ones), then ${\bf A}^{k}{\bf 1}$ converges (after normalization) to the dominant eigenvector of ${\bf A}$. Additionally, $[{\bf A}^{k}{\bf 1}]_{i}$ (i.e., the $i$-th entry of the vector) is equal to the number of length-$k$ paths that emanate from node $i$.

In Fig. 6(b), we plot (as a function of $\omega$) the Pearson correlation coefficient $r$ between node degrees and eigenvector centralities for three cases: (dot-dashed blue curve) intralayer degrees $d_{i}^{(t)}=\sum_{j}A_{ij}^{(t)}$ versus the conditional node centralities $Z_{it}$; (dashed red curve) total degrees $\overline{d}_{i}=\sum_{t}d_{i}^{(t)}$ versus the sum $\sum_{t}Z_{it}(\omega)$ of the conditional node centralities; and (solid gold curve) the degrees $d_{i}^{(2)}=\sum_{j}A_{ij}^{(2)}$ in the Ryanair layer versus $\sum_{t}Z_{it}(\omega)$. As expected, for very small and very large values of $\omega$, the supracentralities correlate strongly with the Ryanair layer and the layer-aggregated network, respectively. Interestingly, for $\omega\approx 0.5$, there is a spike in the correlation between the intralayer degrees and conditional node centralities.

In Sec. SM2 of the Supplementary Materials, we describe the results that we obtain when repeating these computations with PageRank matrices (see Definition 2.6). Figure SM2 gives an interesting contrast to Fig. 6(b). Because PageRank matrices have the same spectral radius, no layer dominates in the limit of small $\omega$, so there is no eigenvector localization (see Thm. 4.5). Instead, in Fig. SM2(b), we observe for small $\omega$ that the conditional centralities correlate strongly with the nodes’ intralayer degrees $d_{i}^{(t)}$. In Sec. SM2, we also compare these findings to results for PageRank “versatility” [22], a generalization of centrality that attempts to quantify important nodes in a multilayer network that may not be particularly important in any individual layer. Because we use an interlayer-adjacency matrix that encodes all-to-all coupling, the different multiplex generalizations of eigenvector centrality that use the Khatri–Rao matrix product [98, 94, 23, 24] are all equal to each other and are also equivalent to the eigenvector centrality that is associated with the layer-averaged adjacency matrix $\sum_{t}\mathbf{A}^{(t)}$. Corollary 4.19 states that this centrality is also equivalent to marginal node centrality for eigenvector supracentralities in the limit $\omega\to\infty$.

5.2 United States Mathematical-Science Program Rankings using a Ph.D. Exchange Network [110]

We apply our supracentrality framework to study the prestige of U.S. mathematical-science doctoral-granting programs by examining a temporal network that encodes the graduation and hiring of Ph.D. recipients in the mathematical sciences. Specifically, we analyze the temporal network from [110], which we constructed using data from the Mathematics Genealogy Project [88]. As in [110], we calculate uniformly and diagonally coupled authority scores, such that a university with a high authority score corresponds to an academic authority. A high-authority university produces desirable students, who tend to be hired by other institutions.

In our study of the Ph.D. exchange network in [110], we restricted our attention to undirected, adjacent-in-time coupling that is encoded by an interlayer-adjacency matrix $\tilde{\bm{A}}$ with entries $\tilde{A}_{tt^{\prime}}=1$ if $|t-t^{\prime}|=1$ and $\tilde{A}_{tt^{\prime}}=0$ otherwise. In the present study, by contrast, we consider the effects of causality by coupling time layers using a directed chain with “layer teleportation”. Specifically, we use an interlayer-adjacency matrix with elements

[TABLE]

Teleportation was introduced for PageRank centrality [33] to allow centrality matrices that are associated with weakly-connected (and even disconnected) networks to satisfy the irreducibility assumptions of Theorems 2.10 and 2.11. Similarly, we use teleportation between layers to satisfy the assumptions of Thm. 3.7; this ensures that the supracentralities are positive and unique. We differentiate between these two types of teleportation by referring to the original PageRank notion as node teleportation and the teleportation in Eq. 29 as layer teleportation. See [111] for an exploration of layer teleportation in further detail.555Theorems 2 and 3 in [111] restate Thms. 4.5 and 4.16 from an earlier draft of the present paper. We have since improved the clarity of these results by stating that $\lambda_{\textrm{max}}$ and $\tilde{\lambda}_{\textrm{max}}$ are “largest positive” eigenvalues.

In Fig. 7, we examine the effect of the layer teleportation parameter $\gamma$ on authority supracentralities. In Figs. 7(a)–(c), we plot the layers’ authority MLCs $x_{t}(\omega)$ (which are given by Definition 3.5) versus the year $t$ for coupling strengths $\omega\in\{10^{1},10^{2},10^{3}\}$ from the weak-coupling, intermediate-coupling, and strong-coupling regimes. See Sec. 3.3 for a description for how we identify coupling regimes. In each panel, we plot the MLCs for three values of the layer teleportation parameter: $\gamma=10^{-2}$, $\gamma=10^{-3}$, and $\gamma=10^{-4}$. In panel (d), we plot $\overline{d}_{t}=\sum_{ij}A_{ij}^{(t)}$, which is the total number of people who earned a Ph.D. in the mathematical sciences in year $t$ who later supervised a graduating Ph.D. student. Observe that $t=1966$ is the year with the largest value of $\overline{d}_{t}$.

In Fig. 7(a) (i.e., for small $\omega$), we observe eigenvector localization onto time layer $t=1982$, whose associated authority matrix has the largest spectral radius among all layers. For $\gamma=10^{-2}$, we also observe a smaller peak at $t=1966$. Comparing Fig. 7(b) to Fig. 7(a) for $\gamma=10^{-2}$, we observe that the peak near $t=1966$ is more pronounced for $\omega=10^{1}$ than for $\omega=10^{2}$. We observed a similar localization phenomenon in [110] for adjacent-in-time coupling. In Fig. 7(c), we illustrate behavior that contrasts starkly with our findings in [110]. Specifically, we see that the interlayer-coupling architecture (which changes with the layer teleportation parameter $\gamma$) has a significant effect on the strong-coupling limit. By varying $\gamma$, one can tune the extent to which older time layers have larger centrality than newer time layers. When $\gamma=10^{-4}$, for example, one can observe in Fig. 7(c) that the MLC of time layers appears to decrease rapidly as $t$ increases. Our asymptotic theory in Sec. 4.2 gives an accurate description of this phenomenon. Observe the dotted black curves, which show the right dominant eigenvector vector $\tilde{\bm{v}}^{(1)}$ of $\tilde{\bm{A}}$; we obtain these MLC curves from Thm. 4.16 in the $\omega\to\infty$ limit.

In Figs. 7(e)–(g), we plot the university rank $r_{i}^{(t)}\in\{1,\dots,N\}$ of Georgia Institute of Technology (GT) that is associated with its conditional centralities for different time layers; we call this its “rank trajectory”. Similar to Figs. 7(a)–(c), these panels show results for $\omega\in\{10,10^{2},10^{3}\}$; in each panel, we again plot results for $\gamma\in\{10^{-2},10^{-3},10^{-4}\}$. In Fig. 7(h), we plot the number of people who earned a Ph.D. in the mathematical sciences from GT who later supervised a graduating Ph.D. student; observe that this increases starting in the 1960s. Starting in the late 1970s, GT’s mathematics program transitioned from being primarily teaching-oriented to being much more research-oriented, with a newly restructured doctoral degree program [26]. We used GT in [110] as a case study to illustrate the methods that we developed in that paper. All centrality trajectories for GT that we present in the present paper differ significantly from those in Taylor et al. [110], who implemented coupled time layers using an undirected chain (i.e., adjacent-in-time coupling with undirected interlayer edges between corresponding universities). In particular, many of the rank trajectories in Figs. 7(e)–(g) suggest that GT has its highest rank in the 1980s, around the time when GT graduated its largest numbers of Ph.D. students who subsequently supervised their own Ph.D. students. We also see that the choices for $\omega$ and $\gamma$ influence the centrality trajectory for GT. We observe that $\gamma$ has a larger effect for the intermediate-coupling and strong-coupling regimes (see panels (f) and (g)) of $\omega$ than it does for the weak-coupling regime. Additionally, for the intermediate choices $\gamma=10^{-3}$ and $\omega=10^{2}$, the rank of GT varies from about 50th to about 15th during when $t\in\{1945,\ldots,1985\}$.

In Sec. SM3 of the Supplementary Materials, we provide additional results for the Ph.D. exchange data. We show in Table SM. 3 that, although $\gamma$ has a strong effect on the MLCs for large values of $\omega$, it does not seem to have a significant effect on top-ranked schools. In Fig. SM3, we show plots that are similar to Figs. 5 and 6, except that the results describe the Ph.D. exchange network, rather than the multiplex airline network. For the limits of small and large $\omega$, we find that the authority supracentralies correlate with the intralayer degrees of a dominating layer ($t=1966$ in this case) and with the nodes’ total degrees, respectively.

6 Conclusions

It is important to develop systematic ways of calculating importances in the form of centralities and their generalizations for nodes, edges, and other structures in multilayer networks. In the present paper, we examined centralities that are based on eigenvectors for two popular classes of multilayer networks: (1) multiplex networks, which encode different types of relationships; and (2) temporal networks, in which relationships change in time. We presented a unifying linear-algebraic framework that generalizes eigenvector-based centralities, such as PageRank and hub/authority scores, for multiplex and temporal networks. A key aspect of our approach involves studying joint, marginal, and conditional centralities that one calculates from the dominant eigenvector of a supracentrality matrix, which couples centrality matrices that are associated with individual layers. See [106] for Matlab code that computes supracentralities and reproduces the experimental results of our paper.

Our main methodological contribution is the extension of the supracentrality framework of [110], which previously was restricted to undirected adjacent-in-time coupling, to more general types of interlayer coupling. Our new, more general framework couples layers through an interlayer-adjacency matrix $\tilde{\bm{A}}$, allowing one to study centralities in multilayer networks with a large family of interlayer coupling topologies. We found that the architecture of $\tilde{\bm{A}}$ significantly impacts supracentralities, and we highlighted that some choices are more appropriate than others for different applications. As an example, in Sec. 5.2, we let $\tilde{\bm{A}}$ encode a directed chain with layer teleportation (see Fig. 1(b) for a visualization) and studied a temporal network that encodes the graduation and hiring of Ph.D. recipients in the mathematical sciences. Our results on this example contrast sharply with those in [110] because of our different choice of interlayer coupling666For an in-depth comparison of the effects directed and undirected coupling between time layers, see our recent book chapter [111]. One of our findings in [111] is that coupling layers using a directed chain (which respects the arrow of time) introduces a bias that increases the centralities of node-layer pairs that are associated with the earliest time layers. For example, we illustrated that one can moderate such a bias using a layer teleportation parameter (and we note that it seems fruitful to study multilayer generalizations of so-called “smart teleportation” [60]).. We also studied a multiplex network that encodes airline transportation in Europe (see Sec. 5.1), where we found a different mechanism that yields boosts in centrality. Specifically, we observed that nodes that are important in both the large-$\omega$ and small-$\omega$ limits can receive centrality boosts for intermediate values of $\omega$ compared to nodes that are important in just one of the limits. For one illustration, see the centrality of the airport LEMB in Table 3 and the black curves in Fig. 5(b).

We explored how different interlayer-coupling architectures (as encoded by $\tilde{\bm{A}}$) and interlayer-coupling strengths (as encoded by $\omega$) influence centralities. Specifically, we identified an interesting interplay between $\tilde{\bm{A}}$ and the architectures of the individual layers in multiplex and temporal networks. To gain insight into this interplay, we performed singular perturbation theory in the limits of weak and strong interlayer coupling (see Sec. 4), which lead to layer decoupling and layer aggregation, respectively. We demonstrated that the limiting supracentralies depend on several factors, including the right and left dominant eigenvectors of $\tilde{\bm{A}}$ and the spectral radii of the layers’ centrality matrices, possibly leading to localization of the dominant eigenvector of a supracentrality matrix onto one or more layers (specifically, the ones whose associated centrality matrices have the largest spectral radii). We focused on studying the eigenvector that is associated with the dominant eigenvalue of a supracentrality matrix, because we obtain the supracentralities from that vector. It may also be interesting to explore other eigenvectors and eigenvalues using the perspective of Jordan decompositions [97].

We expect that our results will be useful not only for centrality analysis, but also for studying matrices that arise in data integration. In the context of our problem, we considered both (i) a set of matrices and (ii) a set of relationships between these matrices. Using a “supramatrix” framework, we constructed and analyzed a matrix that reflects both (i) and (ii). Our perturbative approach for analyzing dominant eigenspaces in the present paper assumes that the matrices are nonnegative and square, but it is not limited to matrices that encode network data. Consequently, we expect our findings to also be insightful in other scenarios that involve combining matrices of data into larger matrices.

One can also use a supracentrality framework to study dynamical processes on multiplex networks. For example, multiplex Markov chains were defined recently [108] using a formulation that was inspired by supracentrality. Such chains can outperform other diffusion models when assessing the importances of nodes in a multiplex network.

Appendix A Proof of Theorem 3.7

To prove the uniqueness and positivity of $\mathbbm{v}(\omega)$, we use the Perron–Frobenius theorem for nonnegative matrices (see Thm. 2.10) [4]. To satisfy the assumptions of Thm. 2.10, we must show that the matrix $\mathbb{C}(\omega)$ is nonnegative and irreducible under our two assumptions: (i) $\tilde{\bm{A}}$ is nonnegative and irreducible; and (ii) the sum $\sum_{t}\mathbf{C}^{(t)}$ is an irreducible nonnegative matrix. By construction, the entries in $\mathbb{C}(\omega)$ are nonnegative, so we only need to prove irreducibility. Because the matrix $\mathbb{C}(\omega)$ is nonnegative, we can interpret it as a weighted adjacency matrix that encodes (possibly) directed and weighted edges between node-layer pairs $\{(i,t)\}$ for $i\in\mathcal{V}=\{1,\dots,N\}$ and $t\in\{1,\dots,T\}$. We will show that the network that is associated with the adjacency matrix $\mathbb{C}(\omega)$ is strongly connected, which implies that it is irreducible.

We start with two observations. First, $\tilde{\bm{A}}$ describes an adjacency matrix for a strongly connected network. Let $\tilde{L}<\infty$ denote the diamater of this network. For any node $k$ and any two layers, $t$ and $t^{\prime}$, it follows that there exists a path from node-layer pair $(k,t)$ to node-layer pair $(k,t^{\prime})$ in the network that is associated with the adjacency matrix $\mathbb{C}(\omega)$. The length of this path is at most $\tilde{L}$. Second, because the matrix $\overline{\mathbf{C}^{(t)}}=T^{-1}\sum_{t}\mathbf{C}^{(t)}$ is irreducible and nonnegative, we can interpret it as an adjacency matrix for a strongly connected network. Let ${L}<\infty$ denote the diameter of this network. For any two nodes, $i$ and $j$, of this network, it follows that there exists a path of length $l\leq L$ from $i$ to $j$. We denote the path by a sequence $\mathcal{P}(i,j)=\{k_{0},k_{1},\dots,k_{l-1},k_{l}\}$ of nodes from $i=k_{0}$ to $j=k_{l}$. We also identify a sequence $\{t_{1},t_{2},\dots,t_{l}\}$ of layers, such that the entry $[\mathbf{C}^{(t_{j})}]_{k_{j-1},k_{j}}$ in matrix $\mathbf{C}^{(t_{j})}$ is positive. For any $j$, there must exist at least one matrix $\mathbf{C}^{(t_{j})}$ for which the $(k_{j-1},k_{j})$-th entry is positive, because the $(k_{j-1},k_{j})$-th entry in $\overline{\mathbf{C}^{(t)}}$ is positive (i.e., because $(k_{j-1},k_{j})$ is an edge) and $\overline{\mathbf{C}^{(t)}}$ is a sum of nonnegative matrices.

For any two node-layer pairs, $(i,s)$ and $(j,t)$, we now prove that there exists a path, with a length of at most $\tilde{L}L$, from $(i,s)$ to $(j,t)$ in the network that is associated with the matrix $\mathbb{C}(\omega)$. We do this by explicitly constructing such a path. We first identify a path $\mathcal{P}(i,j)$ from $i$ to $j$ in the network that is associated with $\overline{\mathbf{C}^{(t)}}$. Consider the following sequence of node-layer pairs:

[TABLE]

where $l\leq{L}$; we define $k_{j}$ and $t_{j}$ as above; and $k_{0}=i$, $t_{0}=s$, $k_{l}=j$, and $t_{l}=t$. By definition, the $(k_{j-1},k_{j})$-th entry in $\mathbf{C}^{(t_{j})}$ is positive for each $j$, implying that the network that is associated with the matrix $\mathbb{C}(\omega)$ has an edge from $(k_{j-1},t_{j})$ to $(k_{j},t_{j})$. We construct a path from $(i,s)$ to $(j,t)$ by taking the sequence in (30) and inserting a path from each term in the sequence to the next term. That is, we insert a path from $(k_{0},t_{0})$ to $(k_{0},t_{1})$ using only node-layer pairs that involve node $k_{0}$. The length of this path is at most $\tilde{L}$. Additionally, from our definition of the path $P(i,j)$, we see that there exists an edge from $(k_{0},t_{1})$ to $(k_{1},t_{1})$. We then insert a path, whose length is also at most $\tilde{L}$, from $(k_{1},t_{1})$ to $(k_{1},t_{2})$ using only node-layer pairs that involve node $k_{1}$. There also exists an edge from $(k_{1},t_{2})$ to $(k_{2},t_{2})$, and so on. We repeat this process until finally we insert a path from $(k_{l-1},t_{l-1})$ to $(k_{l-1},t_{l})$ using only node-layer pairs that involve node $k_{l-1}$, and we note that there exists an edge from $(k_{l-1},t_{l})$ to $(k_{l},t_{l})$. Each of these paths exists because the network that is associated with $\tilde{\bm{A}}$ is strongly connected, and each of these paths has a length of at most $\tilde{L}$.

This construction yields a path from any node-layer pair to any other node-layer pair in the network that is associated with $\mathbb{C}(\omega)$, which proves that the network is strongly connected. Using our construction, we have also obtained an upper bound for the network’s diameter of $\tilde{L}L$. Because $\mathbb{C}(\omega)$ corresponds to a strongly connected network, it is irreducible and nonnegative by the Perron–Frobenius theorem for nonnegative matrices, so the right dominant eigenvector $\mathbbm{v}(\omega)$ is unique and positive. Consequently, the entries $\{W_{it}(\omega)\}$, $\{x_{i}(\omega)\}$, and $\{\hat{x}_{t}(\omega)\}$ are also unique and positive. Because these entries are positive, it follows in turn that $\{Z_{it}(\omega)\}$ and $\{\hat{Z}_{it}(\omega)\}$ are positive and finite.

Finally, if $\mathbb{C}(\omega)$ is also aperiodic, then Thm. 2.10 states that the largest positive eigenvalue of $\mathbb{C}(\omega)$ is larger in magnitude than the other eigenvalues.

Appendix B Proof of Lemma 4.1

We show that each eigenvalue–eigenvector pair of ${\bf C}^{(t)}$ yields an eigenvalue–eigenvector pair of $\mathbb{C}(0)$. Consider the matrix–vector multiplication {{\color[rgb]{0,0,0}\mathbb{C}(0)}}\mathbbm{v}^{(i,t)}, and let k={{\color[rgb]{0,0,0}p}}~{}\text{mod}(N) and t^{\prime}=\lceil{{\color[rgb]{0,0,0}p}}/jN\rceil. We write

[TABLE]

where we use the notation [\hat{\mathbb{C}}\mathbbm{v}^{(i,t)}]_{{\color[rgb]{0,0,0}p}} to denote the {{\color[rgb]{0,0,0}p}}-th entry of the vector {{\color[rgb]{0,0,0}\mathbb{C}(0)}}\mathbbm{v}^{(i,t)}. This implies that {{\color[rgb]{0,0,0}\mathbb{C}(0)}}\mathbbm{v}^{(i,t)}=\mu_{i}^{(t)}\mathbbm{v}^{(i,t)}, so $\mu_{i}^{(t)}$ is an eigenvalue of {{\color[rgb]{0,0,0}\mathbb{C}(0)}} with right eigenvector $\mathbbm{v}^{(i,t)}$. Similarly, {{\color[rgb]{0,0,0}\mathbb{C}(0)}}^{{\color[rgb]{0,0,0}*}}\mathbbm{u}^{(i,t)}=\mu_{i}^{(t)}\mathbbm{u}^{(i,t)}, so $\mathbbm{u}^{(i,t)}$ is the associated left eigenvector.

Appendix C Proof of Theorem 4.5

Equation (17) follows from Lemma 4.1 and Remark 4.4. Let {{\color[rgb]{0,0,0}\mathcal{T}}} be the set of layer indices $\{t\}$ for which $\mu_{1}^{(t)}=\lambda_{\text{max}}(0)$, and define \mathcal{S}=\text{span}\left(\{\mathbbm{v}^{(1,t)}\}_{t\in{{\color[rgb]{0,0,0}\mathcal{T}}}}\right) as the span of the associated eigenvectors (i.e., as the right dominant eigenspace). Because of the continuity of “eigenprojections” (i.e., projections onto eigenspaces) [52] (see Ch. 2), the dominant eigenvector $\mathbbm{v}^{(1)}(\omega)$ converges to the dominant eigenspace $\mathcal{S}$ of $\mathbb{C}(0)$. Similarly, when $\omega\to 0^{+}$, the left dominant eigenvector of $\mathbb{C}(\omega)$ converges to the left dominant eigenspace of $\mathbb{C}(0)$. Given these observations, we only need to

prove that the constants \{\alpha_{{\color[rgb]{0,0,0}t}}\} and \{\beta_{{\color[rgb]{0,0,0}t}}\} satisfy Eqs. (18).

We expand $\lambda_{\textrm{\rm{max}}}(\omega)$, $\mathbbm{u}^{(1)}(\omega)$, and $\mathbbm{v}^{(1)}(\omega)$ for small $\omega$ to obtain order-$k$ approximations:

[TABLE]

We use superscripts to indicate powers of $\omega$ in the terms in the expansion, and we use subscripts for the terms that are multiplied by a power of $\omega$. Note that $\lambda_{0}$, $\mathbbm{v}_{0}$, and $\mathbbm{u}_{0}$, respectively, denote the dominant eigenvalue and its corresponding right and left eigenvectors in the $\omega\to 0^{+}$ limit. Successive terms in the expansions (C) involve higher-order derivatives, and we assume that each term has appropriate smoothness of these functions.

Our strategy is to develop consistent solutions to {\mathbb{C}}(\omega)^{{\color[rgb]{0,0,0}*}}{\mathbbm{u}}^{(1)}(\omega)=\lambda_{\rm{max}}(\omega){\mathbbm{u}}^{(1)}(\omega) and ${\mathbb{C}}(\omega){\mathbbm{v}}^{(1)}(\omega)=\lambda_{\rm{max}}(\omega){\mathbbm{v}}^{(1)}(\omega)$ for progressively larger values of $k$. Let’s consider the equation for the right eigenvector. Starting with the first-order approximation, we insert $\lambda_{\rm{max}}(\omega)\approx\lambda_{0}+\omega\lambda_{1}$ and $\mathbbm{v}^{(1)}(\omega)\approx\mathbbm{v}_{0}+\omega\mathbbm{v}_{1}$ into Eq. (12) and collect the zeroth-order and first-order terms in $\omega$ to obtain

[TABLE]

where $\mathbb{I}$ is the $NT\times NT$ identity matrix. Equation (33) corresponds to the system that is described by Lemma 4.1, implying that the operator $\lambda_{0}\mathbb{I}-\hat{\mathbb{C}}$ is singular and has a {{\color[rgb]{0,0,0}|{{\color[rgb]{0,0,0}\mathcal{T}}}|}}-dimensional null space, where {{\color[rgb]{0,0,0}\mathcal{T}}}=\{t:\mu_{t}=\max_{t}\mu_{t}\} is the set of centrality matrices ${\bm{C}}^{(t)}$ whose largest eigenvalue is equal to the maximum eigenvalue. In particular, $\max_{t}\mu_{t}=\lambda_{0}=\lambda_{\rm{max}}(0)$, and the dominant eigenvectors have a general solution of the form form

[TABLE]

where $\alpha_{t}$ and $\beta_{t}$ are constants that satisfy $1=\sum_{t}\alpha_{t}^{2}=\sum_{t}\beta_{t}^{2}$ to ensure that $\|\bm{\alpha}\|=\|\bm{\beta}\|=1$. (See Lemma 4.1 for definitions of $\mathbbm{u}^{(1,t)}$ and $\mathbbm{v}^{(1,t)}$.)

To determine the vectors \bm{\alpha}=[\alpha_{1},\dots,\alpha_{T}]^{{\color[rgb]{0,0,0}*}} and \bm{\beta}=[\beta_{1},\dots,\beta_{T}]^{{\color[rgb]{0,0,0}*}} of constants that uniquely determine $\mathbbm{u}_{0}$ and $\mathbbm{v}_{0}$, we use Eq. (34) to seek a solvability condition for the first-order terms. We use the fact that the left null space of $\lambda_{0}\mathbb{I}-\hat{\mathbb{C}}$ is the span of $\{\mathbbm{u}^{(1,t)}\}$ to see that [\mathbbm{u}^{(1,t)}]^{{\color[rgb]{0,0,0}*}}\left(\lambda_{0}\mathbb{I}-\hat{\mathbb{C}}\right)\mathbbm{v}_{1}=0 for any $t$. We left-multiply Eq. (34) by [\mathbbm{u}^{(1,t)}]^{{\color[rgb]{0,0,0}*}} and simplify to obtain

[TABLE]

Using the solution of $\mathbbm{v}_{0}$ in Eq. (35), we obtain

[TABLE]

which uses [\mathbbm{u}^{(1,t)}]^{{\color[rgb]{0,0,0}*}}\mathbbm{v}^{(1,t)}=\langle\bm{u}^{(1,t)},\bm{v}^{(1,t^{\prime})}\rangle\delta_{tt^{\prime}}, where $\delta_{ij}$ is the Kronecker delta. We simplify the left-hand side of Eq. (37) to obtain

[TABLE]

The last expression follows from the relations $\hat{\mathbb{A}}=\tilde{\bm{A}}\otimes\bf{I}$, $\mathbbm{v}^{(1,t^{\prime})}={{\bf e}}^{(t^{\prime})}\otimes{\bf v}^{(1,t^{\prime})}$, and $\mathbbm{u}^{(1,t)}={{\bf e}}^{(t)}\otimes{\bf u}^{(1,t)}$, where we recall that ${{\bf e}}^{(t)}$ is a unit vector that consists of zeros in all entries except for entry $t$ (which is a $1$). We equate the expressions (38) and (37) and divide by $\langle{\bf u}^{(1,t)},{\bf v}^{(1,t)}\rangle$ to obtain the equation

[TABLE]

for the right dominant eigenvalue. One proceeds analogously to obtain a similar equation for the left dominant eigenvector, and together these two eigenvector equations yield Eq. (18).

Appendix D Proof of Lemma 4.13

Examining $\tilde{\mathbb{C}}(\epsilon)$, which is given by Eq. (23), yields (using any matrix norm) $\|\tilde{\mathbb{C}}(\epsilon)-\hat{\mathbb{A}}\|=\epsilon\|\hat{\mathbb{C}}\|\to 0^{+}$ as $\epsilon\to 0^{+}$, implying that $\tilde{\mathbb{C}}(0^{+})=\hat{\mathbb{A}}=\tilde{\bm{A}}\otimes\mathbf{I}$. Using the stride permutation $\mathbb{P}$ that is defined by Eq. (24), we write

[TABLE]

Because $\mathbf{I}\otimes\tilde{\bm{A}}$ is block diagonal and each diagonal block is identical, it follows that the spectrum of $\mathbf{I}\otimes\tilde{\bm{A}}$ is identical to that of $\tilde{\bm{A}}$ (although the eigenvalues need to repeat an appropriate number of times), and one can obtain the eigenvectors of the former as functions of the eigenvectors of $\tilde{\bm{A}}$.

Let $\{\tilde{\mu}_{t}\}$ denote the eigenvalues of $\tilde{\bm{A}}$, and let $\tilde{\bm{v}}^{(t)}$ and $\tilde{\bm{u}}^{(t)}$ denote their corresponding right and left eigenvectors, respectively. We now illustrate that $\tilde{\mathbbm{v}}^{(t,j)}=\tilde{{\bf e}}^{(j)}\otimes{\bf v}^{(t)}$ and $\tilde{\mathbbm{u}}^{(t,j)}=\tilde{{\bf e}}^{(j)}\otimes{\bf u}^{(t)}$ are right and left eigenvectors of $\mathbf{I}\otimes\tilde{\bm{A}}$. With {{\color[rgb]{0,0,0}p}}\in\{1,\dots,NT\}, we define t={{\color[rgb]{0,0,0}p}}~{}\text{mod}(T) and k=\lceil t/{{\color[rgb]{0,0,0}p}}T\rceil and obtain

[TABLE]

One can show similarly that (\mathbf{I}\otimes\tilde{\bm{A}})^{{\color[rgb]{0,0,0}*}}\tilde{\mathbbm{u}}^{(t,j)}=\tilde{\mu}_{t}\tilde{\mathbbm{u}}^{(t,j)}, illustrating that $\tilde{\mathbbm{v}}^{(t,j)}$ and $\tilde{\mathbbm{u}}^{(t,j)}$ are right and left eigenvectors that are associated with the eigenvalue $\tilde{\mu}_{t}$ of $\mathbf{I}\otimes\tilde{\bm{A}}$. This implies that \mathbb{P}^{{\color[rgb]{0,0,0}*}}(\bm{A}\otimes\mathbf{I})\mathbb{P}\tilde{\mathbbm{v}}^{(t,j)}=\mu_{t}\tilde{\mathbbm{v}}^{(t,j)}, and left-multiplication by $\mathbb{P}$ gives

[TABLE]

By repeating this procedure using (\tilde{\bm{A}}\otimes\mathbf{I})^{{\color[rgb]{0,0,0}*}} (instead of $\tilde{\bm{A}}\otimes\mathbf{I}$), one can also show that (\tilde{\bm{A}}\otimes\mathbf{I})^{{\color[rgb]{0,0,0}*}}[\mathbb{P}\tilde{\mathbbm{u}}^{(t,j)}]=\tilde{\mu}_{t}[\mathbb{P}\tilde{\mathbbm{u}}^{(t,j)}]. Taken together, this expression and Eq. (D) imply that $\mathbb{P}\tilde{\mathbbm{v}}^{(t,j)}$ and $\mathbb{P}\tilde{\mathbbm{u}}^{(t,j)}$ are right and left eigenvectors of $\hat{\mathbb{A}}=\tilde{\bm{A}}\otimes\mathbf{I}$ that are associated with the eigenvalue $\tilde{\mu}_{t}$. However, for a given value of $t$ (and assuming that the eigenvalues $\{\tilde{\mu}_{t}\}$ are simple), there are $N$ orthogonal right eigenvectors $\{\mathbb{P}\tilde{\mathbbm{v}}^{(t,j)}\}$ and $N$ orthogonal left eigenvectors $\{\mathbb{P}\tilde{\mathbbm{u}}^{(t,j)}\}$ for $j\in\{1,\dots,N\}$. It follows that each eigenvalue $\tilde{\mu}_{t}$ of $\tilde{\mathbb{A}}$ has multiplicity $N$ and associated $N$-dimensional right and left eigenspaces.

Appendix E Proof of Theorem 4.16

We expand $\tilde{\lambda}_{\rm{max}}(\epsilon)$, $\tilde{\mathbbm{v}}^{(1)}(\epsilon)$, and $\tilde{\mathbbm{u}}^{(1)}(\epsilon)$ for small $\epsilon$ to obtain order-$k$ approximations:

[TABLE]

We use superscripts to indicate powers of $\epsilon$ in the terms in the expansion and subscripts for the terms that are multiplied by $\epsilon^{j}$. Note that $\tilde{\lambda}_{0}$, $\tilde{\mathbbm{v}}_{0}$, and $\tilde{\mathbbm{u}}_{0}$, respectively, indicate the dominant eigenvalue and its corresponding right and left dominant eigenvectors in the $\epsilon\to 0^{+}$ limit. Successive terms in the expansions (E) involve higher-order derivatives, and we assume that each term has appropriate smoothness of these functions.

Our strategy is to develop consistent solutions to both of the eigenvalue equations, \tilde{\mathbb{C}}(\epsilon)^{{\color[rgb]{0,0,0}*}}\tilde{\mathbbm{u}}(\epsilon)=\tilde{\lambda}_{\rm{max}}(\epsilon)\tilde{\mathbbm{u}}(\epsilon) and $\tilde{\mathbb{C}}(\epsilon)\tilde{\mathbbm{v}}(\epsilon)=\tilde{\lambda}_{\rm{max}}(\epsilon)\tilde{\mathbbm{v}}(\epsilon)$, for progressively larger values of $k$. Let’s consider the equation for the right eigenvector. Starting with the first-order approximation, we insert $\tilde{\lambda}_{\rm{max}}(\epsilon)\approx\tilde{\lambda}_{0}+\epsilon\tilde{\lambda}_{1}$ and $\tilde{\mathbbm{v}}(\epsilon)\approx\tilde{\mathbbm{v}}_{0}+\epsilon\tilde{\mathbbm{v}}_{1}$ into Eq. (12) and collect the zeroth-order and first-order terms in $\epsilon$ to obtain

[TABLE]

where $\mathbb{I}$ is the $NT\times NT$ identity matrix. Equation (44) corresponds to the system that is described by Lemma 4.1, implying that the operator $\tilde{\lambda}_{0}\mathbb{I}-\hat{\mathbb{A}}$ is singular and has an $N$-dimensional null space. (This is the dominant eigenspace of $\hat{\mathbb{A}}$.) Specifically, Eq. (44) has a general solution of the form

[TABLE]

where $\tilde{\alpha}_{i}$ are constants that satisfy $\sum_{i}\tilde{\alpha}_{i}^{2}=1$ (which implies that $\|\tilde{\mathbbm{v}}_{0}\|_{2}=1$). Additionally,

[TABLE]

where $\sum_{i}\tilde{\beta}_{i}^{2}=1$ (which implies that $\|\tilde{\mathbbm{u}}_{0}\|_{2}=1$).

To determine the vectors \tilde{\bm{\alpha}}=[\tilde{\alpha}_{1},\dots,\tilde{\alpha}_{N}]^{{\color[rgb]{0,0,0}*}} and \tilde{\bm{\beta}}=[\tilde{\beta}_{1},\dots,\tilde{\beta}_{N}]^{{\color[rgb]{0,0,0}*}} of constants that uniquely determine $\tilde{\mathbbm{u}}_{0}$ and $\tilde{\mathbbm{v}}_{0}$, we seek a solvability condition for the first-order terms. Using the fact that the left null space of $\tilde{\lambda}_{0}\mathbb{I}-\hat{\mathbb{A}}$ is the span of $\{\mathbb{P}\tilde{\mathbbm{u}}^{(1,i)}\}$, it follows that [\tilde{\mathbbm{u}}^{(1,i)}]^{{\color[rgb]{0,0,0}*}}\mathbb{P}^{{\color[rgb]{0,0,0}*}}\left(\tilde{\lambda_{0}}\mathbb{I}-\hat{\mathbb{A}}\right)\tilde{\mathbbm{v}}_{1}=0. Therefore, we left-multiply Eq. (45) by [\tilde{\mathbbm{u}}^{(1,i)}]^{{\color[rgb]{0,0,0}*}}\mathbb{P}^{{\color[rgb]{0,0,0}*}} and simplify to obtain

[TABLE]

Using the solution of $\tilde{\mathbbm{v}}_{0}$ in Eq. (E), we obtain

[TABLE]

because \mathbb{P}^{{\color[rgb]{0,0,0}*}}\mathbb{P}=\mathbb{P}\mathbb{P}^{{\color[rgb]{0,0,0}*}}=\mathbb{I} and [\tilde{\mathbbm{u}}^{(1,i)}]^{{\color[rgb]{0,0,0}*}}\tilde{\mathbbm{v}}^{(1,j)}=\langle\tilde{\bm{u}}^{(1)},\tilde{\bm{v}}^{(1)}\rangle\delta_{ij}, where $\delta_{ij}$ is the Kronecker delta. We divide Eq. (49) by $\langle\tilde{\bm{u}}^{(1)},\tilde{\bm{v}}^{(1)}\rangle$ to obtain an $N$-dimensional eigenvalue equation for the dominant eigenvector $\tilde{\bm{\alpha}}$. One can implement the analogous steps for the equations for left dominant eigenvector, and together these two eigenvector equations yield Eq. (27).