Signal Processing Over Multilayer Graphs: Theoretical Foundations and Practical Applications

Signal processing over single-layer graphs has become a mainstream tool owing to its power in revealing obscure underlying structures within data signals. However, many real-life data sets and systems, including those in Internet of Things (IoT), are characterized by more complex interactions among distinct entities, which may represent multilevel interactions that are harder to be captured with a single-layer graph, and can be better characterized by multilayers graph connections. Such multilayer or multilevel data structures can be more naturally modeled by high-dimensional multilayer graphs (MLGs). To generalize traditional graph signal processing (GSP) over MLGs for analyzing multilevel signal features and their interactions, this work proposes a tensor-based framework of MLG signal processing (M-GSP). Specifically, we introduce core concepts of M-GSP and study properties of MLG spectral space, followed by fundamentals of MLG-based filter design. To illustrate novel aspects of M-GSP, we further explore its link with traditional signal processing and GSP. We provide example applications to demonstrate the efficacy and benefits of applying MLGs and M-GSP in practical scenarios.


I. INTRODUCTION
G EOMETRIC signal processing tools have found broad applications in data analysis to uncover obscure or hidden structures from complex datasets [1].Various data sources, such as social networks, Internet of Things (IoT) intelligence, traffic flows and biological images, often feature complex structures that pose challenges to traditional signal processing tools.Recently, graph signal processing (GSP) emerges as an effective tool over the graph signal representation [2].For a signal with N samples, a graph of N nodes can be formed to model their underlying interactions [3].In GSP, a graph Fourier space is also defined from the spectrum space of the representing matrix (adjacency/Laplacian) for signal processing tasks [4], such as denoising [5], resampling [6], and classification [7].Generalization of the more traditional GSP includes signal processing over hypergraphs [8] and simplicial complexes [9], which are suitable to model high-degree multilateral node relationships.
Traditional graph signal processing tools generally describe signals as graph nodes connected by one type of edges.However, real-life systems and datasets may feature multifacet interactions [10].For example, in a video dataset modeled by spatial-temporal graph shown in Fig. 1(a), the nodes may exhibit different types of spatial connections at different temporal steps.It is harder for single-layer graphs to model such multi-facet connections.To model multiple layers of signal connectivity, we explore a high-dimensional graph representation known as multilayer graphs (MLG) [11].
Multilayer graph, also named as multilayer network, is a geometric model containing correlated layers with different structures and physical meanings, unlike traditional singlelayer graphs [10].A typical example is smart grid consisting of two layers shown as Fig. 1(b): the power grid and the computation network.These two layers have different physical connectivity and rules [12].Still, signal interactions across the multiple layers in MLG can be strongly correlated.Thus, separate representations by multiple single layer graphs may fail to capture such characteristics.Consider a network consisting of a physical power layer and a cyber layer, the failure of one layer could trigger the failure of the other [13].One example was the power line damage caused by a storm on September 28th of 2003.Not only did it lead to the failure of several power stations, but also disrupted communications as a result of power station breakdowns that eventually affected 56 million people in Europe [14].
The complexity and multi-level interactions of MLG make the data reside on the irregular and high-dimensional structures, which do not directly lend themselves to standard GSP tools.For example, even though one can represent MLG by a supra-graph unfolding all the layers [15], traditional GSP would treat interlayer and intralayer interactions equivalently in one spectrum space without differentiating the spectra of interlayer and intralayer signal correlations.Recently, there has been growing interest in developing advanced GSP tools to process such multi-level structures.In [16], a joint time-vertex Fourier transform (JFT) is defined to process spatial-temporal graphs by applying graph Fourier transform (GFT) and discrete Fourier transform (DFT).Although JFT can process time-varying datasets, it does not provide more general temporal (interlayer) connectivity in a generic multilayer graph.For flexible interlayer structure, a two-way GFT proposed in [17], defines different graph Fourier spaces for interlayer and intralayer connections, respectively.However, its spatial interactions all lie within a single graph structure, thereby limiting the generalization of MLG.Expanding the works of [17], the tensor-based multi-way graph signal processing framework (MWGSP) of [18] relies on product graph.In MWGSP, separate factor graphs are constructed for each mode of a tensor-represented signal, and a joint spectrum combines factor graph spectra.However, both JFT and MWGSP use the same intralayer connections for all layers.They do not admit different layers with heterogeneous graph structures, thereby limiting the ability to represent a general MLG.
Another challenge in MLG signal processing lies in the need for a suitable mathematical representation.Traditional methods start with connectivity matrices.For example, in [19], a supra-adjacency matrix is defined to represent all layers equivalently while ignoring the natures of different layers.One can also represent each layer with an individual adjacency matrix [20].However, such matrix-based representations mainly focus on the intralayer connections and lack representation for interlayer interactions.A more natural and general way may start with tensor representation [11], which is particularly attractive in handling complex MLG graph analysis.
Our goal is to generalize graph signal processing for multilayer graphs to model, analyze, and process signals based on the intralayer and interlayer signal interactions.To address the aforementioned challenges and to advance MLG processing, we present a novel tensor framework for multilayer graph signal processing (M-GSP).We summarize the main contributions of this work as follows: • Leveraging tensor representation of MLG, we introduce M-GSP from a spatial perspective, in which MLG signals and shifting of signals are defined; • Taking a spectrum perspective, we introduce new concepts of spectrum space and spectrum transform for MLG.For interpretability of spectrum space, we analyze the resulting MLG spectral properties and their distinction from existing GSP tools.• We also present fundamentals of filter design in M-GSP, and suggest several practical applications based on the proposed framework, including those in IoT systems.
We organize the technical presentation as follows.Section II first summarizes preliminaries of traditional GSP and tensor analysis, before presenting representations of multilayer graphs within M-GSP in Section.III.We then introduce the fundamentals of M-GSP spatial and spectrum analysis in Section IV and Section V, respectively.We next discuss MLG filter design in Section VI.We develop the physical insights and spectrum interpretation of M-GSP concepts in Section VII.With the newly proposed M-GSP framework, we provide several example applications to demonstrate its potential in Section VIII, before concluding this paper in Section IX.

II. PRELIMINARIES A. Overview of Graph Signal Processing
Signal processing on graphs [1]- [3] studies signals that are discrete in some dimensions by representing the irregular signal structure using a graph G = {V, F}, where With a graph representation F and a signal vector s, the basic graph filtering (shifting) is defined via The graph spectrum space, also known as the graph Fourier space is defined based on the eigenspace of the representing matrix.Let the eigen-decomposition of F be given by F = VΛV −1 , where V is the matrix with eigenvectors of F as columns, and diagonal matrix Λ consists of the corresponding eigenvalues.The graph Fourier transform (GFT) is defined as whereas the inverse GFT is given by s = Vŝ.
From definitions of GFT, other concepts, such as sampling theory [21], filter design [22], and frequency analysis [4] can be developed for signal processing and data analysis tasks.

B. Introduction of Tensor Basics
Before introducing the fundamentals of M-GSP, we first review some basics on tensors that are useful for multilayer graph analysis.Tensors can be viewed as multi-dimensional arrays.The order of a tensor is the number of indices needed to label a component of that array [23].For example, a thirdorder tensor has three indices.More specially, a scalar is a zeroth-order tensor; a vector is a first-order tensor; a matrix is a second-order tensor; and an M -dimensional array is an M thorder tensor.For convenience, we use bold letters to represent the tensors excluding scalars, i.e., A ∈ R I1×I2•••×I N represents an N th-order tensor with I k being the dimension of the kth order, and use We now start with some useful definitions and tensor operations for the M-GSP framework [23]. 1) Super-diagonal Tensor: 2) Symmetric Tensor: A tensor is super-symmetric if its elements remain constant under index permutation.For example, a third-order In addition, tensors can be partially symmetric in two or more modes as well.For example, a third-order tensor A ∈ R I×I×J is symmetric in the order one and two if A ijk = A jik , for 1 ≤ i, j ≤ I and 1 ≤ k ≤ J.
3) Tensor Outer Product: The tensor outer product between a P th-order tensor U ∈ R The result tensor, whose entries are calculated by The tensor outer product is useful to construct a higher order tensor from several lower order tensors.

4) n-mode Product:
The n-mode product between a tensor U ∈ R I1×•••×I P and a matrix V ∈ R J×In is denoted by Each element in W is given by Note that the n-mode product is a different operation from matrix product.

5) Tensor Contraction:
In M-GSP, the contraction (inner product) between a forth order tensor A ∈ R M ×N ×M ×N and a matrix x ∈ R M ×N in the third and forth order is defined as where y αi = M β=1 N j=1 A αiβj x βj .In addition, the contraction between two fourth-order tensor U, V ∈ R M ×N ×M ×N is defined as whose entries are W αi p = βj U αiβj V βj p .
6) Tensor Decomposition: Tensor decompositions are useful tools to extract the underlying information of tensors.Particularly, CANDECOMP/PARAFAC (CP) decomposition decomposes a tensor as a sum of the tensor outer product of rank-one tensors [23], [24].For example, a third order tensor T ∈ R I×J×K is decomposed by CP decomposition into where integer R denotes the rank, i.e., the lowest number of rank-one tensors in the decomposition.We illustrate CP decomposition for a third-order tensor in Fig. 2(a), which could be viewed as factorization of the tensor.
Another important decomposition is the Tucker decomposition, which is in the form of higher-order PCA.More specifically, Tucker decomposition decomposes a tensor into a core tensor multiplied by a matrix along each mode [23].Defining a core tensor S = [S pqr ] ∈ R P ×Q×R .Defining where The Tucker decomposition for a third-order tensor is illustrated in Fig. 2(b).Tucker decomposition reduces to CP decomposition if the core tensor is limited to be super-diagonal.

III. REPRESENTATIONS OF MLG IN M-GSP
In this section, we introduce the fundamental representations of MLG within the M-GSP framework.

A. Multilayer Graphs
Multilayer graph, also referred to as multilayer network, is an important geometric model in complex systems [10].In this work, we refrain from using the "network" terminology because of its diverse meanings in various field ranging from communication networking to deep learning.From here onward, unless otherwise specified, we shall use the less ambiguous term of multilayer graph.We first provide definitions of multilayer graphs (MLG).
Definition 1 (Multilayer Graph).A multilayer graph with K nodes and M layers is defined as M = {V, L, F}, where denotes the set of layers with each layer being the subsets of V, whereas F is the algebraic representation describing node interactions.
Note that, we mainly focus on the layer-disjoint multilayer graph [10] where each node exists exactly in one layer, since layers denote different phenomena.For example, in a smart grid, a station with functions in both power grid and communication network, is usually modeled as two nodes in a two-layer graph for the network analysis [28].
In multilayer graphs, edges connect nodes in the same layer (intralayer edges) or nodes of different layers (interlayer edges) [29].There are two main types of multilayer graphs: multiplex graph and interconnected graph [30].In a multiplex graph, each layer has the same number of nodes, and each node only connects with their 1-to-1 matching counterparts in other layers to form interlayer connections.Typically, multiplex graphs characterize different types of interactions among the same (or a similar) set of physical entities.For example, the spatial-temporal connections among a set of nodes can be intuitively modeled as a multiplex graph [30].In the interconnected graphs, each layer may have different numbers of nodes without a 1-to-1 counterpart.Their interlayer connections could be more flexible.Examples of a three-layer multiplex graph and a three-layer interconnected graph are shown in Fig. 3, where different colors represent different layers, solid lines represent intralayer connections, and dash lines indicate interlayer connections.

B. Algebraic Representation
To capture the high-dimensional 'multilayer' interactions between different nodes, we use tensor as algebraic representation of MLG for the proposed M-GSP framework [11].
1) MLG with same number of nodes in each layer: To better interpret the tensor representation of a multilayer graph, we start from a simpler type of MLG, in which each layer contains the same number of nodes.For a multilayer graph M = {V, L} with |L| = M layers and N nodes in each layer, i.e., |l i | = N for 1 ≤ i ≤ M , it could be interpreted as embedding the interactions between a set of N 'entities' (not nodes) into a set of M layers.The nodes in different layers can be viewed as the projections of the entities.For example, the video datasets could be modeled by the spatial connections between objects (entities) into different temporal frames (layers).
Mathematically, the process of embedding (projecting) entities can be viewed as a tensor product, and graph connections can be represented by tensors [11].For convenience, we use Greek letters α, β, • • • to indicate each layer and Latin letters i, j, • • • to indicate each interpretable 'entity' with corresponding node in each layer.Given a set of entities whose sole nonzero element is its i−th element (equal to 1) to characterize each entity i.Thus, interactions of two entities can be represented by a second-order tensor A X = N i,j=1 a ij e i • e j ∈ R N ×N , where a ij is the intensity of the relationship between entity i and j.Similarly, given a set of layers L = {l 1 , l 2 , • • • , l M }, a vector e α ∈ R M can capture the properties of the layer α, and the connectivity between two layers could be represented by Following this approach, connectivity between the projected nodes of the entities in the layers can be represented by a fourth-order tensor (11) where w αiβj is the weight of connection between the entity i's projected node on layer α and the entity j's projected node on layer β.More specially, the fourth-order tensor becomes the adjacency tensor of the multilayer graph, where each entry A αiβj = w αiβj characterizes the edge between the entity i's projected node on layer α and the entity j's projected node on layer β.Thus, similar to the adjacency matrix whose 2-D entries indicate whether and how two nodes are pairwise connected by a simple edge in the normal graphs, we adopt an adjacency tensor A ∈ R M ×N ×M ×N to represent the multilayer graph with the same number of nodes in each layer as follows.
Definition 2 (Adjacency Tensor).A multilayer graph M, with |L| = M layers and |l i | = N nodes in each layer i, can be represented by a fourth-order adjacency tensor A ∈ R M ×N ×M ×N defined as Here, each entry A αiβj of the adjacency tensor A indicates the intensity of the edge between the entity j's projected node on layer β and entity i's projected node on layer α.
Clearly, for a single layer graph, e α is a scalar 1 and the fourth-order tensor degenerates to the adjacency matrix of a normal graph.Similar to A ij in an adjacency matrix which indicates the direction from the node v j to v i , A αiβj also indicates the direction from the node v βj to the node v αi in a MLG.Note that, vectors e i and e α are not eigenvectors of the adjacency tensor.They are merely the vectors characterizing features of the entities and layers, respectively.We shall discuss the MLG-based spectrum space in Section V.
Given an adjacency tensor, we can define the Laplacian tensor of the multilayer graph similar to that in a singlelayer graph.Denoting the degree (or multi-strength) of the entity's i's projected node v αi on layer α as d(v αi ) which is a summation over weights of different natures (inter-and intra-layer edges), the degree tensor D ∈ R M ×N ×M ×N is defined as a diagonal tensor with entries D αiαi = d(v αi ) for 1 ≤ i ≤ N, 1 ≤ α ≤ M , whereas its other entries are zero.The Laplacian tensor can be defined as follows.The Laplacian tensor can be useful to analyze propagation processes such as diffusion or random walk [11].Both adjacency and Laplacian tensors are important algebraic representations of the MLG depending on datasets and user objectives.For convenience, we use a tensor F ∈ R M ×N ×M ×N as a general representation of a given MLG M. As the adjacency tensor is more general, the representing tensor F refers to the adjacency tensor in this paper unless specified otherwise.
2) Representation of General MLG: Representing a general multilayer graph with different number of nodes in each layer always remains a challenge if one aims to distinguish the interlayer and intralayer connection features.In JFT [16] and MWGSP [18], all layers must reside on the same underlying graph structure which restrict the number of nodes to be the same in each layer.Similarly, a reconstruction is also needed to represent a general MLG by the forth-order tensor in M-GSP.Note that, although M-GSP also needs a reconstruction to represent a general MLG, we allow different layers with heterogeneous graph structures, which provides additional flexibility than JFT and MWGSP.
There are mainly two ways to reconstruct: 1) Add isolated nodes to layers with fewer nodes to reach N nodes [12] and set the augmented signals as zero; and 2) Aggregate several nodes into super-nodes for layers with |l i | > N [31] and merge the corresponding signals.Since isolated nodes do not interact with any other nodes, it does not change the topological structure of the original multilayer architecture in the sense of signal shifting while the corresponding spectrum space could still be changed.The aggregation method depends on how efficiently we can aggregate redundant or similar nodes.Different methods can be applied depending on specific tasks.For example, if one wants to explore the cascading failure in a physical system, the method based on isolated nodes is more suitable.For the applications, such as video analysis where pixels can be intuitively merged as superpixels, the aggregation method can be also practical.
In addition, although the fourth-order representing tensor can be viewed as the projection of several entities into different layers in Eq. ( 11), the entities and layers can be virtual and not necessarily physical to capture the underlying structures of the datasets.The information within the multilayer graphs, together with definitions of the underlying virtual entities and layers, should only depend on the structure of the multilayer graphs.We will illustrate this further in Section VII-B.

C. Flattening and Analysis
In this part, we introduce the flattening of the multilayer graph, which could simplify some operations in the tensorbased M-GSP.For a multilayer graph M = {V, L, F} with M layers and N nodes in each layer, its fourth-order representing tensor F ∈ R M ×N ×M ×N can be flattened into a second-order matrix to capture the overall edge weights.There are two main flattening schemes in the sense of entities and layers, respectively: • Layer-wise Flattening: The representing tensor F can be flattened into • Entity-wise Flattening: The representing tensor F can be flattened into These two flattening methods provide two ways to interpret the graph structure.In the first method, the flattened multilayer graph has M clusters with N nodes in each cluster.The nodes in the same cluster have the same function (belong to the same layer).In the second method, the flattened graph has N clusters with M nodes in each cluster.Here, the nodes in the same cluster are from the same entity.Examples of the tensor flattening of a two-layer graph with 3 nodes in each layer are shown in Fig. 4. From the examples, we can see that the diagonal blocks in R N ×N are the intralayer connections for each layer and other blocks describe the interlayer connections through layer-wise flattening; and the diagonal block in R M ×M describe the 'intra-entity' connections and other elements represent the 'inter-entity' connections in entity-wise flattening.Although these two flattening schemes define the same MLG with a different indexing of vertices, they are still helpful to analyze the MLG spectrum space.

IV. SPATIAL DEFINITIONS IN M-GSP
Based on the tensor representation, we now define signals and signal shifting over the multilayer graphs.In GSP, each signal sample is the attribute of one node.Typically, a graph signal can be represented by an N -length vector for a graph with N nodes.Recall that in traditional GSP [3], basic signal shifting is defined with the representing matrix as the shifting filter.Thus, in M-GSP, we can also define the signals and signal shifting based on the filter implementation.
In M-GSP, each signal sample is also related to one node in the multilayer graph.Intuitively, if there are K = M N nodes, there are M N signal samples in total.Similar to GSP, we use the representing (adjacency/Laplacian) tensor F ∈ R M ×N ×M ×N as the basic MLG-filter.Since the input signal and the output signal of the MLG-filter should be consistent in the tensor size, we define a special form of M-GSP signals to work with the representing tensor as follows.
where the entry s αi is the signal sample in the projected node of entity i on layer α.
Note that, if the multilayer graph degenerates to a singlelayer graph with M = 1, the multilayer graph signal becomes an N -length vector, which is similar to that in the traditional GSP.Similar to the representing tensor, the tensor signal s ∈ R M ×N can also be flattened as a vector in R M N : • Layer-wise flattening: Given the definitions of multilayer graph signals and filters, we now introduce the definitions of signal shifting in M-GSP.In traditional GSP, the signal shifting is defined as product between signal vectors and representing matrix.Similarly, we define the shifting in the multilayer graph based on the contraction (inner product) between the representing tensor and tensor signals.
Definition 5 (Signal Shifting over Multilayer Graphs).Given the representing tensor F ∈ R M ×N ×M ×N and the tensor signal s ∈ R M ×N defined over a multilayer graph M, the signal shifting is defined as the contraction (inner product) between F and s in one entity-related order and one layerrelated order, i.e., where is the contraction between F and s defined in Eq. (7).
The elements in the shifted signal s are calculated as From Eq. ( 17), there are two important factors to construct the shifted signal: 1) The signal in the neighbors (both intra-and inter-layer interactions) of the node v αi ; and 2) The intensity of interactions between the node v αi and its neighbors.Then, the signal shifting is related to the diffusion process over the multilayer graphs.More specifically, if F is the adjacency tensor, signals shift in directions of edges.To better illustrate the signal shifting based on the adjacency tensor, we use a two-layer directed network shown in Fig. 5 as an example.In this multilayer graph, the original signal s is defined as and the adjacency tensor A is defined as for each fiber.Then, the shifted signal is calculated as From Eq. ( 19), we can see that the signal shift one step following the direction of the links.Meanwhile, if F is the Laplacian tensor, Eq. ( 17) can be written as which is the weighted average of difference with neighbors.

V. MULTILAYER GRAPH SPECTRUM SPACE
In traditional GSP, graph spectrum space is defined according to the eigenspace of the representing matrix [3].Similarly, we define the MLG spectrum space based on the decomposition of the representing tensor.Since tensor decomposition is less stable when exploring the factorization of a specific order or when extracting the separate features in the asymmetric tensors, we will mainly focus on spectral properties of undirected multilayer graphs in this section for simplicity and clarity of presentation.For directed MLG, we provide alternative spectrum definitions in the Appendix and will explore detailed analysis in future works.
A. Joint Spectral Analysis in M-GSP For a multilayer graph M = {V, L, F} with M layers and N nodes, the eigen-tensor V ∈ R M ×N of the representing tensor F is defined in the tensor-based multilayer graph theory [11] as More specifically, F ∈ R M ×N ×M ×N can be decomposed as where λ k is the eigenvalues and V k ∈ R M ×N is the corresponding eigen-tensor.Note that V αi just relabels the index of V k , and there is no specific order for V αi here.Similar to the traditional GSP where the graph Fourier space is defined by the eigenvectors of the representing matrix, we define the joint MLG Fourier space as follows.
Definition 6 (Joint Multilayer Graph Fourier Space).For a multilayer graph M = {V, L, F} with M layers and N nodes, the MLG Fourier space is defined as the space consisting of all spectral tensors {V 1 , • • • , V M N }, which characterizes the joint features of entities and layers.
Recall that in GSP, the GFT is defined based on the inner product of V −1 and the signals s defined in Eq. ( 2).Similarly, we can define the M-GFT based on the spectral tensors of the representing tensor F to capture joint features of inter-and intra-layer interactions as follows.
The joint M-GFT (M-JGFT) can be defined as the contraction between U F and the tensor signal s ∈ R M ×N , i.e., Now, we show how to obtain the eigen-tensors.Implementing the flattening analysis, we have the following properties.
Property 1.The two types of flattened tensor in Eq. ( 13) and Eq. ( 14) lead to the same eigenvalues.
Proof.Suppose (λ, x) is an eigenpair of A F L , i.e., we have Thus, λ is also an eigenvalue of A F N .
This property shows that the flattened tensors are the reshaped original representing tensor, and could capture some of the spectral properties as follows.
Property 2. Given the eigenpair (λ F L , x) of the layerwise flattened tensors, the eigenpair (λ, V) of the original representing tensor can be calculated as λ = λ F L , and V αi = x N (α−1)+i .Similarly, given the eigenpair (λ F N , y) of the entity-wise flattened tensor, the eigenpair (λ, V) of the original representing tensor can be calculated as λ = λ F N , and V αi = y M (i−1)+α .
The Property 2 shows that we can calculate the eigentensor from the flattened tensor to simplify the decomposition operations.Moreover, the M-JGFT is the bijection of GFT in the flattened MLG, with vertices indexed by both the layers and the entities.However, such M-JGFT analyzes the inter-and intra-layer connections jointly while ignoring the individual features of entities and layers.Next, we will show how to implement the order-wise frequency analysis in M-GSP based on tensor decomposition.

B. Order-wise Spectral Analysis in M-GSP
In an undirected multilayer graph, the representing tensor (adjacency/Laplacian) F is partially symmetric between orders one and three, and between orders two and four, respectively.Then, the representing tensor can be written with the consideration of the multilayer graph structure under orthogonal CP-decomposition [26] as follows: where The CP decomposition factorizes a tensor into a sum of component rank-one tensors, which describe the underlying features of each order.Although approximated algorithms are implemented to obtain the optimal decomposition, CP decomposition still achieves great success in real scenarios, such as feature extraction [32] and tensor-based PCA analysis [33].A detailed discussion of tensor decomposition and its implementation in M-GSP are provided in Section VII-D.In Eq. ( 29), f α and e i capture the features of layers and entities, respectively, which can be interpreted as the subspaces of the MLG.More discussions about the frequency interpretation of order-wise M-GSP spectrum and connections to MWGSP spectrum are presented in Section VII-C.
Note that, if there is only one layer in the multilayer graph, Eq. ( 29) reduces to the eigen-decomposition of a normal single-layer graph, i.e., F = N i=1 λ i e i • e i With the decomposed representing tensor in Eq. ( 29), the order-wise MLG spectrum is defined as follows.
Definition 8 (Order-wise MLG Spectral Pair).For a multilayer graph M = {V, L, F} with M layers and N nodes, the order-wise MLG spectral pairs are defined by {λ αi , f α , e i }, With the definition of order-wise MLG spectral pair, we now explore their properties.Considering Ṽαi = f α • e i , we have the following property, which indicates the availability of a joint MLG analysis based on order-wise spectrum.
Property 3. The factor tensor Ṽαi of the representing tensor F is the approximated eigen-tensor of F.
Proof.Suppose that Ṽαi is one factor tensor of F obtained from Eq. (29).
Let δ[k] denote the Kronecker delta.Since f α forms an orthonormal basis, then the inner product would satisfy Similarly, Then, we have Thus, which indicates Then, Ṽαi is the approximated eigen-tensor.
This property indicates the relationship between the orderwise MLG spectral pair and the joint eigen-tensors.
This property generalizes the orthogonality of the spectral tensor into a similar definition of matrix.
We now introduce the order-wise MLG spectral transform.Similar to Eq. ( 24), the joint transform can be defined as Note that each element of ŝ in Eq. ( 37) can be calculated as Let with each element Clearly, the M-GFT can be obtained as Then, we have the following definition of M-GFT based on order-wise spectrum.
Definition 9 (Order-wise M-GFT).Given the spectral vectors the layer-wise M-GFT can be defined as and the entity-wise M-GFT can be defined as The general M-GFT based on order-wise spectrum is defined as If there is only one layer in the multilayer graph, the M-GFT calculated with s T ∈ R N as (ŝ N ) T = (sE e ) T ∈ R N , which has the same form as the traditional GFT in Eq. ( 2).
In addition, since f α and e i are orthonormal basis of undirected MLG, the inverse M-GFT can be calculated as Different from joint MLG Fourier space in Section V-A, the order-wise MLG spectrum provides an individual analysis on layers and entities separately, and a reliable approximated analysis on the underlying MLG structures jointly.To minimize confusion, we abbreviate joint M-GFT in Eq. ( 24) as M-JGFT.M-GFT refers to the order-wise transform in Eq. ( 46) in the remaining parts if there is no specification.

C. MLG Singular Tensor Analysis
In addition to the eigen-decomposition, the singular value decomposition (SVD) is another important decomposition to factorize a matrix.In this part, we provide the higher-order SVD (HOSVD) [25] of the representing tensor as an alternative definition of spectrum for the multilayer graphs.
Given the multilayer graph M = {V, L, F} with M layers and N nodes in each layer, its representing tensor F ∈ R M ×N ×M ×N can be decomposed via HOSVD as where In ] is a unitary (I n × I n ) matrix, with I 1 = I 3 = M and I 2 = I 4 = N .S is a complex (I 1 × I 2 × I 3 × I 4 )-tensor of which the subtensor S in obtained by fixing nth index to α have is the n-mode singular value, and U (i) are the corresponding n-mode singular vectors.For an undirected multilayer graph, the representing tensor is symmetric for every 2-D combination.Thus, there are two modes of singular spectrum, i.e., (γ α , f α ) for mode 1, 3, and (σ i , e i ) for mode 2, 4.More specifically, U (1) = U (3) = (f α ) and U (2) = U (4) = (e i ).Since the joint singular tensor captures the consistent information of entities and layers, it can be calculated as Note that the diagonal entries of S are not the eigenvalues or frequency coefficients of the representing tensor in general.
The multilayer graph singular space is defined as follows.
Definition 10 (Multilayer Graph Singular Space).For a multilayer graph M = {V, L, F} with M layers and N nodes, the MLG singular space is defined as the space consisting of all singular tensors { V1 Similar to order-wise spectral analysis in Section V-B, we can define the MLG singular tensor transform (M-GST) based on the singular tensors as follows.
Definition 11 (M-GST).Suppose that U s = (f α • e i ) ∈ R M ×N ×M ×N consists of the singular vectors of the representing tensor F in Eq. (48), where The M-GST can be defined as the contraction between U s and the tensor signal s ∈ R M ×N , i.e., š = U s s. (50) If the singular vectors are included in , the layer-wise M-GST can be defined as and the entity-wise M-GST can be defined as Inverse M-GST can be defined similarly as in Eq. ( 47) with unitary W e and W f .Compared to the eigen-tensors in Eq. ( 22), the singular tensors come from the combinations of the singular vectors, thus are capable of capturing information of layers and entities more efficiently.Eigen-decomposition, however, focuses more on the joint information and approximate the separate information of layers and entities.We shall provide further discussion on the performance of different decomposition methods in Section VII-D.The intuition of applying HOSVD in MLG analysis and its correlations to GSP are also provided in Section VII-A3.

D. Spectrum Ranking in the Multilayer Graph
In traditional GSP, the frequencies are defined by the eigenvalues of the shift, whereas the total variation is an alternative measurement of the order of the graph frequencies [3].Similarly, we use the total variation of λ αi based on the spectral tensors to rank the MLG frequencies.Let |λ| max be the joint eigenvalue of the adjacency tensor A with the largest magnitude.The M-GSP total variation is defined as follows: where || • || 1 is the element-wise l 1 norm.Other norms could also be used to define the total variation.For example, the l 2 norm could be efficient in signal denoising [3].The graph frequency related to λ αi is said to be a higher frequency if its total variation T V (V αi ) is larger, and its corresponding spectral tensor V αi is a higher frequency spectrum.If the representation tensor refers to Laplacian tensor, i.e., L = D − A, the frequency order is in contract to the adjacency tensor A as GSP [3].We shall provide more details on interpretation of MLG frequency in Section VII-A.

VI. FILTER DESIGN
In this section, we introduce an MLG filter design together with its properties based on signal shifting.

A. Polynomial Filter Design
Polynomial filters are basic filters in GSP [7], [34].In M-GSP, first-order filtering consists of basic signal filtering, i.e., Similarly, a second order filter can be defined as additional filtering on first-order filtered signal, i.e., whose entries s αi are calculated as where is the contraction defined in Eq. ( 8).
Property 5.The τ -th order basic shifting filter f τ (s) can be calculated as This property is the M-GSP counterpart to the traditional linear system interpretation that complex exponential signals are eigenfunctions of linear systems [3], and provide a quicker implementation of higher-order shifting.With the k-order polynomial term, the adaptive polynomial filter is defined as where {α k } are parameters to be estimated from data.Adaptive polynomial filter is useful in semi-supervised classification [35] and exploits underlying geometric topologies.We will illustrate further and provide application examples based on MLG polynomial filtering in Section VIII.

B. Spectral Filter Design
Filtering in the graph spectrum space is useful in GSP frequency analysis.For example, ordering the Laplacian graph spectrum in a descent order by the graph total variation [3], i.e., high frequency to low frequency, the GFT of s ∈ R N is calculated as ŝ = V T G s.By removing k elements in the low frequency part, i.e., ŝ = [ŝ 1 , • • • , ŝN−k , 0, • • • , 0], a high-pass filter can be designed as where Similarly, in M-GSP, a spectral filter is designed by filtering in the spectrum space together with the inverse M-GFT.With Eq. ( 46) and Eq. ( 50), spectral filtering of s is defined as where functions g(•) and f (•) are designed by the specific tasks.For example, if one wants to design a layerwise filter capturing the smoothness of signals in the MLG singular space, the function g(•) can be designed as g by ordering the layer-wise singular vectors in the descent order of singular values.

C. Discussion
We briefly discuss the interpretation of polynomial filters and spectral filters.From the spatial perspective, MLG polynomial filter is a weighted sum of messaging passing on the multilayer graph in different orders, shown as Eq.(66).Each node collects information from both inter-and intra-layer neighbors, before combining them with its own information.From the spectrum perspective, M-GSP polynomial filters are eigenfunctions of linear systems, which are special cases of M-GSP spectral filters shown in Eq. (69) The M-GSP spectral filters assign different weights to each M-GSP spectrum via functions f (•) and g(•) depending on specific tasks.Both M-GSP polynomial filters and spectral filters can be useful for high-dimensional IoT signal processing.More discussions and examples of M-GSP filters are presented in Section VIII.

VII. DISCUSSION AND INTERPRETATIVE INSIGHTS A. Interpretation of M-GSP Frequency 1) Interpretation of Graph Frequency:
To better understand its physical meaning, we start with the total variation in digital signal processing (DSP).The total variation in DSP is defined as differences among signals over time [36].Moreover, the total variations of frequency components have a 1-to-1 correspondence to frequencies in the order of their values.If the total variation of a frequency component is larger, the corresponding frequency with the same index is higher.This means that, a higher frequency component changes faster over time and exhibits a larger total variation.Interested readers could refer to [3], [8] for a detailed interpretation of total variation in DSP.Now, let us elaborate the graph frequency motivated by the cyclic graph.Rewrite the finite signals in DSP as vectors, i.e., the signal shifting can be interpreted as the shift filtering corresponding to a cyclic graph shown in Fig. 6.Suppose that its adjacency matrix is written as Then, the shifted signal over the cyclic graph is calculated as T , which shifts the signal at each node to its next node.More specifically, C N can be decomposed as C N = VΛV −1 where the eigenvalues is the discrete Fourier matrix.Inspired by the DSP, the eigenvectors in V are the spectral components (spectrum) of the cyclic graph and the eigenvalues are related to the graph frequencies [3].
Generalizing the adjacency matrix of the cyclic graph to the representing matrix F M of an arbitrary graph, the graph Fourier space consists of the eigenvectors of F M and the graph frequencies are related to the eigenvalues.More specifically, the graph Fourier space can be interpreted as the manifold or spectrum space of the representing matrix.As aforementioned, the total variations of frequency components reflect the order of frequencies, we can also use the total variation, i.e., T V (e i ) = ||e i − 1 |λ|max e i || 1 , to rank the graph frequencies, where e i is the spectral component related to the eigenvalue λ i in F M .Similar to DSP, the graph frequency indicates the oscillations over the vertex set, i.e., how fast the signals change over the graph shifting.
2) Interpretation of MLG Frequency: Now, return to M-GSP.Given spectral tensors V k ∈ R M ×N of a multilayer graph, a signal s ∈ R M ×N can be written in a weighted sum of the spectrum, i.e., s = k a k V k .Viewing the spectral tensor as a signal component, the total variation is in the form of differences between the original signal and its shifted version in Eq. ( 53).If the signal component changes faster over the multilayer graph, the corresponding total variation is larger.Since we relate higher frequency component with a larger total variation, the MLG frequency indicates how fast the signal propagates over the multilayer graph under the representing tensor.If a signal s contains more high frequency components, it changes faster under the representing tensor.
3) Interpretation of MLG Singular Tensors: As discussed in Section VII-A1, the name of graph Fourier space arises from the adjacency matrix of the cyclic graph.However, when the algebra representation is generalized to an arbitrary graph, especially the Laplacian matrix, the definition of graph spectrum is less related to the Fourier space in DSP but can be interpreted as the manifold or subspace of the representing matrix instead.In literature, SVD is an efficient method to obtain the spectrum for signal analysis, such as spectral clustering [37] and PCA analysis [38].It is straightforward to generalize graph spectral analysis to graph singular space, especially for the Laplacian matrix.In M-GSP, the order-wise singular vectors can be interpreted as subspaces characterizing features of layers and entities, respectively.Since HOSVD is robust and efficient, transforming signals to the MLG singular space (M-GST) for the analysis of underlying structures can be a useful alternative for M-GFT.

B. Interpretation of Entities and Layers
To gain better physical insight of entities and layers, we discuss two categories of datasets: • In most of the physical systems and datasets, signals can be modeled with a specific physical meaning in terms of layers and entities.In smart grid, for example, each station can be an entity, connected in two layers of computation and power transmission, respectively.Another example is video in which each geometric pixel point is an entity and each video frame form a layer.Each layer node denotes the pixel value in that video frame.M-GSP can be intuitive tool for these datasets and systems.• In some scenarios, however, the datasets usually only has a definition of layers without meaningful entities.In particular, for multilayer graphs with different numbers of nodes, we may insert some isolated artificial nodes to augment the multilayer graph.Often in such applications, it may be harder to identify the physical meaning of entities.Here, the entities may be virtual and are embedded in the underlying structure of the multilayer graph.
Although definition of a virtual entity may vary with the chosen adjacency tensor, it relates to the topological structure in terms of global spectral information.For example, in Fig. 7, we can use two different definitions of virtual entities.Although the representing tensors for these two definitions differ, their eigenvalues remain the same.Considering also layer-wise flattening, the two supra-matrices are related by reshaping, by exchanging the fourth and fifth columns and rows.They still have the same eigenvalues, whereas the eigentensors can also be the same by implementing the reshaping operations.Note that, to capture distinct information from entities, their spectra would change with different definitions of virtual entities.
C. Distinctions from Existing GSP Works 1) Graph Signal Processing: Generally, M-GSP extends traditional GSP into multilayer graphs.Although one can stack all MLG layers to represent them with a supra-matrix, such matrix representation makes GSP inefficient in extracting features of layers and entities separately.Given a supra-matrix of the MLG, the layers of nodes can not be identified directly from its index since all the nodes are treated equivalently.However, the tensor representation provides clear identifications on layers in its index.Moreover, in GSP, we can only implement a joint transform to process inter-and intralayer connections together, while the M-GSP provide a more flexible choice on joint and order-wise analysis.In Section V-A, the joint M-GSP analysis introduced can be viewed as the bijection of GFT in the flattened MLG, with vertices indexed by both layers and entities.Beyond that, we flexibly provide order-wise spectral analysis based on tensor decompositions, which allow the order-wise analysis on layers and nodes.One can select suitable MLG-based tools depending on tasks.The joint spectral analysis can be implemented if we aim to explore layers and entities fully, whereas the order-wise spectral and singular analysis are more efficient in characterizing layers and entities separately.
2) Joint Time-Vertex Fourier Analysis: In [16], a joint timevertex Fourier transform (JFT) is defined by implementing GFT and DFT consecutively.As discussed in Section VII-A1, the time sequences can be interpreted under a cyclic graph, and thus reside on a MLG structure.However, JFT assumes that all the layers have the same intra-layer connections, which limits the generalization of the MLG analysis.Differently, the tensorbased representation allows heterogeneous structures for the intra-layer connections, which makes M-GSP more general.
3) Multi-way Graph Signal Processing: In [18], MWGSP has been proposed to process high-dimensional signals.Given Kth-order high-dimensional signals, one can decompose the tensor signal in different orders, and construct one graph for each.Graph signal is said to reside on a high-dimensional product graph obtained by the product of all individual factor graphs.Although the MW-GFT is similar to M-GFT for K = 2, there still are notable differences in terms of spectrum.First, MWGSP can only process signals without exploiting a given structure since multiple graph spectra would arise from each order of the signals.For a multilayer graph with a given structure, such as physical networks with heterogeneous intralayer connections, MWGSP does not naturally process it efficiently and cohesively.The order-wise spectra come from factor graphs of each order in MWGSP while M-GSP spectra are calculated from the tensor representation of the whole MLG.Second, MWGSP assumes all the layers residing on a homogeneous factor graph and restricts the types of manageable MLG structure.For example, in a spatial temporal dataset, a product graph, formed by the product of spatial connections and temporal connections, assumes the same topology in each layer.However, many practical systems and datasets feature more complex geometric interactions.M-GSP provide a more intuitive and natural framework for such MLG.In summary, despite some shared similarities between MW-GFT and M-GFT in some scenarios, they serve different purposes and are suitable for different underlying data structures.

D. Comparison of Different Decomposition Methods
To compare recovery accuracy of representing tensor using different tensor decomposition methods, we examine eigentensor decomposition, HOSVD, optimal CP decomposition and Tucker decomposition in MLGs randomly generated from the Erdös−R ėnyi (ER) random graph ER(p, q, M, N ).Here M is the number of layers with N nodes in each layer, p is the intralayer connection probability and q is the interlayer connection probability.We apply different decomposition methods of similar complexity, and compute errors between decomposed and original tensors.From Table I, we can see that the eigen-tensor decomposition and HOSVD exhibit better overall accuracy.Generally, eigen-tensor decomposition is better suited for applications emphasizing joint features of layers and entities.On the other hand, HOSVD is effective at separating individual features of layers and entities.Note that, in addition to recovery accuracy, different decompositions may have different performance advantages when capturing different data features that can be better measured with different metrics.

VIII. APPLICATION EXAMPLES
We now provide some illustrative application examples within our M-GSP framework.

A. Analysis of Cascading Failure in IoT Systems
Analysis of cascading failure in IoT systems based on the spreading processes in multilayer graphs has recently attracted significant interests [39].Modeling the failure propagation in complex systems by shifting over multilayer graph, M-GSP spectrum theory can help the analysis of system stability.In this part, we introduce a M-GSP analysis for cascading failure over multilayer cyber-physical systems based on epidemic model [12].Shown in Fig. 1(b), a cyber-physical system with M layers and N nodes in each layer can be intuitively modeled by a MLG with adjacency tensor A ∈ R M ×N ×M ×N .
Here, we consider the susceptible-infectious-susceptible (SIS) model [40] for the failure propagation.In the SIS model, each node has two possible states: susceptible (not fail) or infectious (fail).At each time slot, the infectious node may cause failure to other nodes through directed links at certain infection rates, or it may heal itself spontaneously at a selfcure rate.The initial attack make several nodes infectious.
Since the nodes in the same layer correspond to the same functionality, e.g., power transmission, nodes on the same layer have the same self-cure rate and infection rate.The notations of the epidemic model are given as follows: • θ αβ : infection rate describing failure propagation probability from nodes on layer β to those on layer α; • P i,α (t): failure probability of the projected node of entity i on layer α at time t; • i,α (t): transition probability that the projected node of entity i on layer α shifts from infectious state to susceptible state at time t; • σ i,α (t): transition probability that the projected node of entity i on layer α remains susceptible at time t.Since an infectious node becomes susceptible if it cures itself without being infected by its neighbors, we have Similarly, a susceptible node remains susceptible without being infected by its neighbors.Thus, The state transition forms a Markov chain, for which we derive the failure probability as where We can define a transition tensor T ∈ R M ×N ×M ×N with elements T βj αi to characterize the failure propagation in Eq. (73).In steady state, P i,α (τ ) = P i,α (τ − 1).Moreover, where Pi,α is the failure probability of the projected node of entity i on layer α in steady state.Following [12], we can arrive that Pi,α = 0 if the spectral radius of the transition tensor ρ(T) < 1, which indicates no failed nodes in steady state.Thus, ρ(T) could serve as an indicator for system robustness.Here, to avoid being repetitive, we merely introduce a simple example of MLG-based cascading failure analysis.Interested readers may refer to our work [12] for more details.With a better understanding of M-GSP spatial shifting, one can develop more general analysis for various failure models in the multilayer IoT systems.

B. Spectral Clustering
Clustering is a widely used tool in a variety of applications such as social network analysis, computer vision, and IoT.Spectral clustering is popular and effective among many variants.Modeling dataset by a normal graph before spectral clustering, significant performance improvement is possible in structured data [37].In this part, we introduce M-GSP spectral clustering and demonstrate its application in RGB image segmentation.
To model an RGB image using MLG, we can directly treat its three colors as three layers.To reduce the number of nodes for computational efficiency, we first build N superpixels for a given image and represent each superpixel as an entity in the 7: Assign all pixels inside one superpixel to the cluster of that superpixel; 8: Output: The segmented image.multilayer graph.Here, we define the feature of a superpixel according to its RGB pixel values.For interlayer connections, each node connects with its counterparts in other layers.For intralayer connections on layer , we calculate the Gaussianbased distance between two superpixels according to if |s i ( ) − s j ( )| 2 ≤ t ; otherwise, W ij = 0, where s i ( ) is the superpixel value on layer , δ is an adjustable parameter and t is a predefined threshold.Different layers may have different structures.The threshold t is set to be the mean of all pairwise distances.As such, an RGB image is modeled as multiplex graph with M = 3 and N nodes.We now consider MLG-based spectral clustering.For image segmentation, we focus on the properties of entities (i.e., superpixels), and implement spectral clustering over entitywise spectrum by proposing Algorithm 1.The previous discussions have been summarized in steps 1-3.In Step 4, different schemes may be used to calculate spectrum, including spectral vector via tensor factorization in Eq. ( 29), and singular vector in Eq. (48).Step 5 determines K based on the largest arithmetic gap in eigenvalues.Traditional clustering methods, such as k-means clustering [37], can be carried out in Step 6.
To test the proposed Algorithm 1, we first compare its results with those from GSP-based method and traditional kmeans clustering by using a public BCCD blood cell dataset shown in Fig. 8(a).In this dataset, there are mainly three types of objects, i.e., White Blood Cell (WBC) vs. Red Blood Cell (RBC) vs. Platelet (P).We set the number of clusters to 3 and N = 1000.For GSP-based spectral clustering, we construct graphs based on the Gaussian model by using information from all 3 color values   In addition to visual inspection of results for such images, we are further interested to numerically evaluate the performance of the proposed methods against some state-of-art methods for several more complex datasets that contain more classes.For this purpose, we test our methods on the BSD300 and BSD500 datasets [41].BSD300 contains 300 images with labels, and BSD500 contains 500 images with labels.We first cluster each image, and label each cluster with the best map of cluster orders against the ground truth.Numerically, we use mIOU (mean of Intersection-over-Union), also known as the mean Jaccard Distance, for all clusters in each image to measure the performance.The Jaccard Distance between two groups A and B is defined as A larger mIOU indicates stronger performance.To better illustrate the results, we considered two setups of datasets, i.e., one containing fewer classes (coarse) and one containing all images (all).We compare our methods together with k-means, GSP-based spectral clustering, invariant information clustering (IIC) [42], graph-based segmentation (GS) [43], back propagation (BP) [44] and differentiable feature clustering (DFC) [45].The best performance is marked in bold.From the results of Table II

C. Semi-Supervised Classification
Semi-supervised classification is an important practical application for IoT intelligence.In this application, we apply MLG polynomial filters for semi-supervised classification.Traditional GSP defines adaptive filter as where W is an adjacency matrix based on pairwise distance or a representing matrix constructed from the adjacency matrix.
where M (•) is a mapping of filtered signals to their discrete labels.For example, in a {±1} binary classification, one can assign a label to a filtered signal against a threshold (e.g.0).Some other objective functions include labeling uncertainty, Laplacian regularization, and total variation.Using estimated parameters, we can filter the signal one more time to determine labels for some unlabeled data by following the same process.
Similarly, in an MLG, we can also apply polynomial filters for label estimation.Given an arbitrary dataset X = [x 1 , • • • , x N ] ∈ R K×N with N signal points and K features for each node, we can construct a MLG by defining M = K layers based on features and N entities based on signal points.The inter-and intra-layer connections are calculated by the Gaussian distance with different parameters.Let its adjacency tensor A ∈ R M ×N ×M ×N .A signal is defined by which is an extended version of graph signal.Note that we do not necessarily need to order signals by placing zeros in the rear.We only write the signal as Eq. ( 79) for notational convenience.We now apply polynomial filters on signals, i.e., and For a filtered signal s X ∈ R M ×N (X = 1, 2), we define a function to map 2-D signals into 1-D by calculated the columnwise mean of s X , i.e., Next, we can define a function M (•) on sX and consider certain objective functions in filter design.To validate the efficacy of polynomial filtering in the MLG framework, we test h 1 (•) and h 2 (•) for the binary classification problem on the Cleveland Heart Disease Dataset.In this dataset, there are  297 data points with 13 feature dimensions.We directly build a MLG with N = 297 nodes in each of the M = 13 layers.More specifically, we directly use the labels as s.For h 1 (•) (AF), we set a i = 0 for at least one i > 0. Using MSE as objective function, we apply a greedy algorithm to estimate parameters {a i }.We limit the highest polynomial order to 10.For h 2 (•) (APF), we estimate a classification threshold via the mean of sX by setting the polynomial order i = 10.
We compare our methods with GSP-based method in similar setups as in aforementioned examples.The only difference is that we use sX in M-GSP and use s = f ([s T L 0 T UL ] T ) in GSP for mapping and classification.We also present the results of label propagation and SVM for comparison.We randomly split the test and training data for 100 rounds.From the results shown in Fig. 9, GSP-based and M-GSP based methods exhibit better performance than traditional learning algorithms, particularly when the fraction of test samples is large.In general, M-GSP based methods demonstrate superior performance among all methods owing to its strength to extract 'multilayer' features, which could potentially benefit semisupervised classification tasks in IoT systems.

D. Dynamic Point Cloud Analysis
3D perception plays an important role in the high growth fields of IoT devices and cyber-physical systems, and continues to drive many progresses made in advanced point cloud processing [46].Here, we propose a short time M-GST method to analyze dynamic point cloud.Given a dynamic point cloud with M frames and at most N points in each frame, we model it as a multilayer graph with M layers and N nodes in each layer.More specifically, we test the singular spectrum analysis over the motion sequences of subject 86 in the CMU database [47].To implement the M-GST, we first divide the motion sequence into several shorter sequences with N f frames.Next for each shorter sequence, we model interlayer connections by connecting points with the same label among successive frames.For points in the same frame, we connect two points based on the Gaussian-kernel within a Euclidean threshold τ s [6].Let x i be the 3D coordinates of the ith point.We assign an edge weight between two points x i and x j as a nonzero A ij = exp(− x i −x j 2 2 /σ 2 ) only if x i −x j 2 2 ≤ τ s .Next, we estimate the spatial and temporal basis vectors of each shorterterm sequences by HOSVD in Eq. (48).Finally, we use the 3D coordinates of all points in each shorter-term sequences as signals and calculate their M-GST.To illustrate the results of M-GST, we examine the spectrogram similar to that of shorttime Fourier transform (STFT) [48].In Fig. 10, we transform the signal defined by the coordinates in Z dimension via M-GST and illustrate the transformation results for the divided frame sequence.From Fig. 10, one can easily identify different motions based on the MLG singular analysis.
To explore motions in dynamic point clouds, we can also apply the entity-wise MLG highpass filters described in Section VI-B to capture critical details of human bodies.More specifically, we select the first 140 frames in 'walking' and define the norm of three coordinates as signals.We select 5 body joints (entities) in each temporal frame (layers) shown as Fig. 11.From the results shown, entity 1 and entity 2 exhibit periodic patterns which are linked to the leg motion.Entity 3 (head) shows little movement relative to the main body.Entity 4 and entity 5 (hands) display more irregular patterns since they do not directly identify 'walking'.To summarize, the MLG highpass filter can efficiently capture some key information of body movement and identify the meaning of nodes (entities).These and related information can assist further analysis of dynamic point clouds including compression and classification.
Our future works shall target more practical applications of point cloud on IoT devices, including point cloud compression, low-complex point cloud segmentation and robust denoising.

E. Other Potential Applications in IoT Systems
Along with the widespread deployment of IoT technologies, system structures become increasingly complex.Traditional graph-based tools are less adept at modeling 'multilayer' graph interactions.The more general model of M-GSP provides additional opportunities for IoT applications.Here, we suggest several potential scenarios in IoT systems for M-GSP: • IoT networks with multilayer structure fit naturally to MLG, for which M-GSP can be designed for various tasks such as intrusion detection, resource management and state prediction; • Because of the dynamic nature in practical IoT networking, even signals on single-layer IoT systems naturally fit a spatial-temporal graph model, which can be also characterized by MLG.For such dynamic IoT systems, M-GSP tools, such as adaptive filters and MLG learning machines, can be developed for signal prediction and node classification.
Overall, the power of MLG in extracting underlying 'multilayer/multi-level' structures in the IoT systems makes M-GSP a potentially important tool in handling highdimensional signal processing and learning tasks.

IX. CONCLUSION
In this work, we present a novel tensor-based framework of multilayer graph signal processing (M-GSP) that naturally generalizes the traditional GSP to multilayer graphs.We first present the basic foundation and definitions of M-GSP including MLG signals, signal shifting, spectrum space, singular space, and filter design.We also provide interpretable discussion and physical insights through numerical results and examples to illustrate the strengths, general insights, and benefits of novel M-GSP framework.We further demonstrate exciting potentials of M-GSP in data processing applications through experimental results in several practical scenarios.
With recent advances in tensor algebra and multilayer graph theory, more opportunities are emerging to explore M-GSP and its applications.One such interesting problem is the efficient construction of multilayer graph, where M-GSP spectrum properties could improve the robustness of graph structure.Another promising direction is to develop multilayer graph neural networks based on the M-GSP spectral convolution.Additional future directions include the development of M-GSP sampling theory and fast M-GFT.

APPENDIX
Unlike for undirected graphs, representing tensors of directed graphs is asymmetric, thereby making each layer or entity characterized by a pair of spectral vectors.To find the spectrum space of a directed multilayer graph, we also present two ways to compute: • Flattening analysis: Similar to the representing tensor of undirected graphs, we flatten the representing tensor as a second-order supra-matrix, and define spectrum space as the reshaped eigenvectors of the supra-matrix.The flattened matrix A F X (or A F N , A F L ) can be decomposed as where E ∈ R M N ×M N is the matrix of eigenvectors and Σ = diag(λ i ) is a diagonal matrix of eigenvalues.Then, we can reshape the eigenvectors, i.e., each column of E as V k ∈ R M ×N , and reshape each row of E −1 as U k ∈ R M ×N .Consequently, the original tensor can be decomposed into • Tensor Factorization: We can also compute the spectrum from the tensor factorization based on CP-decomposition where R is the rank of tensor, a k , c k ∈ R M characterize layers, b k , d k ∈ R N characterize entities, and characterize the joint features.Since there are M N nodes, it is clear that R ≤ M N .Note that, for a single layer, Eq. (85) reduces to Moreover, if V = (v k ) and U = (u T k ) = V −1 are orthogonal bases, Eq. ( 87) is in a consistent form of the eigen-decomposition in a single-layer normal graph.In addition, Eq. ( 28) is also a special case of Eq. (85) if the multilayer graph is undirected.Since tensor decomposition is less stable when exploring the factorization of a specific order or when extracting the separate features in the asymmetric tensors, we will defer more general analysis of directed networks to future works.

Fig. 1 .
Fig. 1.Multilayer Graphs and Applications: (a) Video: each layer represents one frame of the video and the edges capture the spatial-temporal relationships; (b) Cyber-Physical System (CPS): each layer represents one component in CPS and the edges capture the physical connections.
and F ∈ R N ×N is the representing matrix (e.g., adjacency/Laplacian) describing the geometric structure of the graph G. Graph signals are the attributes of nodes that underlie the graph structure.A graph signal can be written as vector s = [s 1 , s 2 , • • • , s N ] T ∈ R N where the superscript (•) T denotes matrix/vector transpose.

Definition 3 (
Laplacian Tensor).A multilayer graph M, with |L| = M layers and |l i | = N nodes in each layer i, can be represented by a fourth-order Laplacian tensor L ∈ R M ×N ×M ×N defined as L = D − A, where A is the adjacency tensor and D is the degree tensor.

Definition 4 (
Signals over Multilayer Graphs).For a multilayer graph M = {V, L, F}, with |L| = M layers and |l i | = N nodes in each layer i, the definition of multilayer graph signals is a second-order tensor
characterize features of layers and entities, respectively.

Algorithm 1
MLG-based Unsupervised Image Segmentation 1: Input: RGB Image I ∈ R P ×Q×3 ; 2: Build N superpixels for the image I and calculate the value of superpixel based on the mean of all pixels inside that superpixel, i.e., s ∈ R N ×3 ; 3: Construct a multilayer graph A ∈ R M ×N ×M ×N ; 4: Find entity-wise spectrum E = [e 1 , • • • , e N ] ∈ R N ×N ; 5: Select the first K important leading spectrum based on the eigenvalues (singular values) of E as C ∈ R N ×K ; 6: Cluster each row of C, and assign the ith superpixel into jth cluster if the ith row of C is clustered into jth group;

3 =1
|s i ( ) − s j ( )| 2 to form edge connections in a single graph.There is only a single δ and t in the Gaussian model.For M-GSP, we use the MLG singular vectors (MLG-SVD), and tensor factorization (MLG-FZ) for spectral clustering, separately.Their respective results are compared in Fig 8. WBCs are marked yellow, and RBCs are marked green.Platelet (P) is marked blue.From the illustrative results, MLG methods exhibit better robustness and are better in detecting regions under noise.Comparing

Fig. 8 .
Fig. 8. Example of BCCD Datasets and Segmented Images: (a) the original image; (b)-(e) segmented images under different methods (WBCs are marked yellow, RBCs are marked green, and Platelet (P) is marked blue).

Fig. 10 .
Fig. 10.Example of Transformed Signals in a Dynamic Point Cloud.

TABLE I ERROR
OF DECOMPOSING THE REPRESENTING TENSOR

TABLE II RESULTS
OF IMAGE SEGMENTATION IN IMAGE BSD300 AND BSD500 [45] can see that larger number of clusters in the first two columns generate worse performance.There are two natural reasons.First, the mapping of the best order of cluster labels is more difficult for more classes.Second, the graph-based spectral clustering is sensitive to the number of K leading spectra and the structure of graphs.Regardless, MLGbased methods still demonstrate better performance.Moreover, even when we use the same total number of nodes in a single layer graph and multilayer graph for another fairness comparison in terms of complexity, i.e., N = 300 for graph and N = 100 for MLG, MLG-based methods still perform better than graph-based methods in this example application.MLG methods have competitive performances compared to the state-of-the-art methods.Note that, under proper training, neural network (NN)-based methods may give good results in cases with many clusters as suggested in[45].