# Introducing Graph Signal Processing over Multilayer Networks: Theoretical
Foundations and Frequency Analysis

Songyang Zhang, Qinwen Deng, and Zhi Ding,  S. Zhang, Q. Deng and Z. Ding are
with Department of Electrical and Computer Engineering, University of
California, Davis, CA, 95616. (E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>and zding@ucdavis.edu).

###### Abstract

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 datasets and systems are characterized by more complex
interactions among distinct entities, which may represent multi-level
interactions that are harder to be modeled with a single-layer graph, and can
be better captured by multiple layers of graph connections. Such
multilayer/multi-level data structure can be more naturally modeled with high-
dimensional multilayer network (MLN). To generalize traditional graph signal
processing (GSP) over multilayer networks for the analysis of multi-level
signal features and their interactions, we propose a tensor-based framework of
multilayer network signal processing (M-GSP) in this work. Specifically, we
introduce the core concepts of M-GSP and study the properties of MLN spectral
space. We then present the fundamentals of MLN-based filter design. To
illustrate the novelty of M-GSP, we further describe its connections to
traditional digital signal processing and GSP. We provide example applications
to demonstrate the efficacy and benefits of applying multilayer networks and
the M-GSP in practical scenarios.

###### Index Terms:

Multilayer network, graph signal processing, tensor, data analysis

##  I Introduction

Geometric 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, 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
spectral 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 multi-lateral 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 multi-facet interactions [10]. For example, in a video dataset
modeled by spatial temporal graph shown in Fig. 1, 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 networks [11].

Multilayer network (MLN) is a geometric model containing correlated layers
with different structures and physical meanings, unlike traditional single-
layer graphs. A typical example is smart grid consisting of two layers shown
as Fig. 1: 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 MLN 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].

Figure 1: Multilayer Networks 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.

The complexity and multi-level interactions of MLN 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
MLN by a supra-graph unfolding all the layers [15], traditional GSP would
treat interlayer and intralayer interactions equivalently in one spectral
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 two-step
transform is proposed to process spatial-temporal graphs. Graph Fourier
transform (GFT) is applied first in the spatial domain (intralayer) and then
in the temporal domain (interlayer). In this framework, different graph
Fourier spaces are defined for interlayer and intralayer connections
respectively. However, all the spatial interactions reside within a single
graph structure which limits the generalization of MLN. In [17], a joint time-
vertex Fourier transform (JFT) is defined by implementing GFT and DFT
consecutively. Although JFT can process the time-varying datasets, it does not
process more general temporal (interlayer) connectivity in a generic
multilayer network. Alternatively, a tensor-based multi-way graph signal
processing framework (MWGSP) relies on product graph [18]. In this framework,
separate factor graphs are constructed for each mode of a tensor-represented
signal, and a joint spectrum combines the spectra of all factor graphs. Since
MWGSP focuses on product of all factor graphs estimated from signals, it is
not well suited to a multilayer network with a given structure.

Another challenge in MLN 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
[20]. A more natural and general way may start with tensor representation
[11], which is particularly attractive in handling complex MLN graph analysis.

Our goal is to generalize graph signal processing for multilayer networks to
model, analyze, and process signals based on the intralayer and interlayer
signal interactions. To address the aforementioned challenges and to advance
MLN processing, we present a novel tensor framework for graph signal
processing over multilayer networks (M-GSP). We summarize the main
contributions of this work as follows:

  * •

Leveraging tensor representation of MLN, we introduce the definitions of
signals and shifting over MLN.

  * •

We define new concepts of spectral space and spectrum transform for M-GSP. For
interpretability of the spectral space, we analyze the resulting MLN spectral
properties and their distinctions from existing GSP tools.

  * •

We also present fundamentals of filter design in M-GSP, and suggest several
practical applications based on the proposed framework.

We organize the technical contents as follows. Section II first summarizes the
preliminaries of traditional GSP and tensor analysis. We then introduce the
fundamentals of M-GSP and frequency analysis in Section III and Section IV,
respectively. We next present MLN filter design in Section V. We provide the
physical insights and spectrum interpretation of M-GSP concepts in Section VI.
Within the proposed M-GSP framework, we present several example applications
and demonstrate its effectiveness in Section VII, before summarizing our
conclusions in Section VIII.

##  II Preliminaries

###  II-A Overview of Graph Signal Processing

Signal processing on graphs [1, 3, 2] studies signals that are discrete in
some dimensions by representing the irregular signal structure using a graph
$\mathcal{G}=\\{\mathcal{V},\mathbf{F}\\}$, where
$\mathcal{V}=\\{v_{1},v_{2},\cdots,v_{N}\\}$ is a set of $N$ nodes, and
$\mathbf{F}\in\mathbb{R}^{N\times N}$ is the representing matrix (e.g.,
adjacency/Laplacian) describing the geometric structure of the graph
$\mathcal{G}$. Graph signals are the attributes of nodes that underlie the
graph structure. A graph signal can be written as vector
$\mathbf{s}=[s_{1},s_{2},\cdots,s_{N}]^{\mathrm{T}}\in\mathbb{R}^{N}$ where
the superscript $()^{\mathrm{T}}$ denotes matrix/vector transpose.

With a graph representation $\mathbf{F}$ and a signal vector $\mathbf{s}$, the
basic graph filtering is defined via $\mathbf{s}^{\prime}=\mathbf{Fs}$.

The graph spectral space, also known as the graph Fourier space is defined
based on the eigenspace of the representing matrix. Let the eigen-
decomposition of $\mathbf{F}$ be given by
$\mathbf{F}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1}$, where $\mathbf{V}$ is
the matrix with eigenvectors of $\mathbf{F}$ as columns, and diagonal matrix
$\mathbf{\Lambda}$ consists of the corresponding eigenvalues. The graph
Fourier transform (GFT) is defined as

$\hat{\mathbf{s}}=\mathbf{V}^{-1}\mathbf{s},$ (1)

whereas the inverse GFT is given by $\mathbf{s}=\mathbf{V}\hat{\mathbf{s}}$.

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.

###  II-B Introduction of Tensor Basics

Before introducing the fundamentals of M-GSP, we first review some basics on
tensors that are useful for multilayer network 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 third-order
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$th-order tensor. For convenience, we use bold
letters to represent the tensors excluding scalars, i.e.,
$\mathbf{A}\in\mathbb{R}^{I_{1}\times I_{2}\cdots\times I_{N}}$ represents an
$N$th-order tensor with $I_{k}$ being the dimension of the $k$th order, and
use $A_{{i_{1}}\cdots i_{N}}$ to represent the entry of $\mathbf{A}$ at
position $(i_{1},i_{2},\cdots,i_{N})$ with $1\leq i_{k}\leq I_{k}$ in this
work. If $\mathbf{A}_{f}$ has a subscript $f$, we use $[A_{f}]_{{i_{1}}\cdots
i_{N}}$ to denote its entries.

We now start with some useful definitions and tensor operations for the M-GSP
framework [23].

####  II-B1 Super-diagonal Tensor

An $N$th-order tensor $\mathbf{A}\in\mathbb{R}^{I_{1}\times I_{2}\cdots\times
I_{N}}$ is super-diagonal if its entries $A_{i_{1}i_{2}\cdots i_{N}}\neq 0$
only for $i_{1}=i_{2}=\cdots=i_{N}$.

####  II-B2 Symmetric Tensor

A tensor is super-symmetric if its elements remain constant under index
permutation. For example, a third-order $\mathbf{A}\in\mathbb{R}^{I\times
I\times I}$ is super-symmetric if
$A_{ijk}=A_{jik}=A_{kij}=A_{kji}=A_{jik}=A_{jki}$. In addition, tensors can be
partially symmetric in two or more modes as well. For example, a third-order
tensor $\mathbf{A}\in\mathbb{R}^{I\times I\times J}$ is symmetric in the order
one and two if $A_{ijk}=A_{jik}$, for $1\leq i,j\leq I$ and $1\leq k\leq J$.

####  II-B3 Tensor Outer Product

The tensor outer product between a $P$th-order tensor
$\mathbf{U}\in\mathbb{R}^{I_{1}\times\cdots\times I_{P}}$ with entries
$U_{i_{1}...i_{P}}$ and a $Q$th-order tensor
$\mathbf{V}\in\mathbb{R}^{J_{1}\times\cdots\times J_{Q}}$ with entries
$V_{j_{1}...j_{Q}}$ is denoted by

$\mathbf{W}=\mathbf{U}\circ\mathbf{V}.$ (2)

The result $\mathbf{W}\in\mathbb{R}^{I_{1}\times\cdots\times I_{P}\times
J_{1}\times\cdots\times J_{Q}}$ is a $(P+Q)$th-order tensor, whose entries are
calculated by $W_{i_{1}...i_{P}j_{1}...j_{Q}}=U_{i_{1}...i_{P}}\cdot
V_{j_{1}...j_{Q}}.$ The tensor outer product is useful to construct a higher
order tensor from several lower order tensors.

####  II-B4 n-mode Product

The n-mode product between a tensor
$\mathbf{U}\in\mathbb{R}^{I_{1}\times\cdots\times I_{P}}$ and a matrix
$\mathbf{V}\in\mathbb{R}^{J\times I_{n}}$ is denoted by

$\mathbf{W}=\mathbf{U}\times_{n}\mathbf{V}\in\mathbb{R}^{I_{1}\times\cdots\times
I_{n-1}\times J\times I_{n+1}\times\cdots\times I_{P}}.$ (3)

Each element in $\mathbf{W}$ is given by $W_{i_{1}i_{2}\cdots
i_{n-1}ji_{n+1}\cdots i_{P}}=\sum_{i_{n}=1}^{I_{n}}U_{i_{1}\cdots
i_{P}}V_{ji_{n}}.$ Note that the $n$-mode product is a different operation
from matrix product.

####  II-B5 Tensor Contraction

In M-GSP, the contraction (inner product) between a forth order tensor
$\mathbf{A}\in\mathbb{R}^{M\times N\times M\times N}$ and a matrix
$\mathbf{x}\in\mathbb{R}^{M\times N}$ in the third and forth order is defined
as

$\mathbf{y}=\mathbf{A}\diamond\mathbf{x}\in\mathbb{R}^{M\times N},$ (4)

where $y_{\alpha i}=\sum_{\beta=1}^{M}\sum_{j=1}^{N}A_{\alpha i\beta
j}x_{\beta j}$.

In addition, the contraction between two fourth-order tensor
$\mathbf{U},\mathbf{V}\in\mathbb{R}^{M\times N\times M\times N}$ is defined as

$\mathbf{W}=\mathbf{U}\odot\mathbf{V}\in\mathbb{R}^{M\times N\times M\times
N},$ (5)

whose entries are $W_{\alpha i\epsilon p}=\sum_{\beta j}U_{\alpha i\beta
j}V_{\beta j\epsilon p}$.

####  II-B6 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].
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]. Other
typical decompositions include Higher-Order SVD (HOSVD) [25], orthogonal CP-
decomposition [26], and Tensor-Train decomposition [27]. Interested readers
are referred to the tutorial [23] for more details. Additional examples and
illustrations of tensor decomposition in M-GSP are also provided in the
Appendix B because of page limit.

##  III Fundamentals of M-GSP

In this section, we introduce the basic definitions in the proposed M-GSP
framework.

Figure 2: Example of multilayer networks: (a) A three-layer interconnected
network; (b) a three-layer multiplex network.

###  III-A Multilayer Network

Before introducing the foundations of M-GSP, we first provide definitions of
multilayer networks (MLN) [23].

######  Definition 1 (Multilayer Network).

A multilayer network with $K$ nodes and $M$ layers is defined as
$\mathcal{M}=\\{\mathcal{V},\mathcal{L},\mathbf{F}\\}$, where
$\mathcal{V}=\\{v_{1},v_{2},\cdots,v_{K}\\}$ is the set of nodes,
$\mathcal{L}=\\{l_{1},l_{2},\cdots,l_{M}\\}$ denotes the set of layers with
each layer $l_{i}=\\{v_{i_{1}},\cdots,v_{i_{n}}\\}$ being the subsets of
$\mathcal{V}$, whereas $\mathbf{F}$ is the algebraic representation describing
node interactions.

Note that, we mainly focus on the layer-disjoint multilayer network [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 network for the network analysis [28].

In multilayer networks, edges connect nodes in the same layer (intralayer
edges) or nodes of different layers (interlayer edges) [29]. There are two
main types of multilayer networks: multiplex network and interconnected
network [30]. In a multiplex network, 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 networks
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 network [30]. In the
interconnected networks, 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 network and a three-layer
interconnected network are shown in Fig. 2, where different colors represent
different layers, solid lines represent intralayer connections, and dash lines
indicate interlayer connections.

###  III-B Algebraic Representation

To capture the high-dimensional ‘multilayer’ interactions between different
nodes, we use tensor as algebraic representation of MLN for the proposed M-GSP
framework [10].

####  III-B1 MLN with same number of nodes in each layer

To better interpret the tensor representation of a multilayer network, we
start from a simpler type of MLN, in which each layer contains the same number
of nodes. For a multilayer network $\mathcal{M}=\\{\mathcal{V},\mathcal{L}\\}$
with $|\mathcal{L}|=M$ layers and $N$ nodes in each layer, i.e., $|l_{i}|=N$
for $1\leq i\leq 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 network connections can be represented by tensors
[11]. For convenience, we use Greek letters $\alpha,\beta,\cdots$ to indicate
each layer and Latin letters $i,j,\cdots$ to indicate each interpretable
‘entity’ with corresponding node in each layer. Given a set of entities
$\mathcal{X}=\\{x_{1},x_{2},\cdots,x_{N}\\}$, one can construct a vector
$\mathbf{e}_{i}\in\mathbb{R}^{N}$ to characterize each entity $i$. Thus,
interactions of two entities can be represented by a second-order tensor
$\mathbf{A}_{X}=\sum_{i,j=1}^{N}a_{ij}\mathbf{e}_{i}\circ\mathbf{e}_{j}\in\mathbb{R}^{N\times
N}$, where $a_{ij}$ is the intensity of the relationship between entity $i$
and $j$. Similarly, given a set of layers
$\mathcal{L}=\\{l_{1},l_{2},\cdots,l_{M}\\}$, a vector
$\mathbf{e}_{\alpha}\in\mathbb{R}^{M}$ can capture the properties of the layer
$\alpha$, and the connectivity between two layers could be represented by
$\mathbf{A}_{L}=\sum_{\alpha,\beta=1}^{M}b_{\alpha\beta}\mathbf{e}_{\alpha}\circ\mathbf{e}_{\beta}\in\mathbb{R}^{M\times
M}$. Following this approach, connectivity between the projected nodes of the
entities in the layers can be represented by a fourth-order tensor

$\mathbf{A}=\sum_{\alpha,\beta=1}^{M}\sum_{i,j=1}^{N}w_{\alpha i\beta
j}\mathbf{e}_{\alpha}\circ\mathbf{e}_{i}\circ\mathbf{e}_{\beta}\circ\mathbf{e}_{j}\in\mathbb{R}^{M\times
N\times M\times N},$ (6)

where $w_{\alpha i\beta j}$ is the weight of connection between the entity
$i$’s projected node in layer $\alpha$ and the entity $j$’s projected node in
layer $\beta$. More specially, if we select the vector
$\mathbf{e}_{i}=[0,\cdots,0,1,0,\cdots,0]^{\mathrm{T}}$ in which the only
nonzero element is the $i$th element (equal to 1) for both layers and
entities, the fourth-order tensor becomes the adjacency tensor of the
multilayer network, where each entry $A_{\alpha i\beta j}=w_{\alpha i\beta j}$
characterizes the edge between the entity $i$’s projected node in layer
$\alpha$ and the entity $j$’s projected node in layer $\beta$. 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 $\mathbf{A}\in\mathbb{R}^{M\times N\times M\times N}$ to
represent the multilayer network with the same number of nodes in each layer
as follows.

######  Definition 2 (Adjacency Tensor).

A multilayer network $\mathcal{M}$, with $|\mathcal{L}|=M$ layers and
$|l_{i}|=N$ nodes in each layer $i$, can be represented by a fourth-order
adjacency tensor $\mathbf{A}\in\mathbb{R}^{M\times N\times M\times N}$ defined
as

$\mathbf{A}=(A_{\alpha i\beta j}),\quad 1\leq\alpha,\beta\leq M,1\leq i,j\leq
N.$ (7)

Here, each entry $A_{\alpha i\beta j}$ of the adjacency tensor $\mathbf{A}$
indicates the intensity of the edge between the entity $j$’s projected node in
layer $\beta$ and entity $i$’s projected node in layer $\alpha$.

Clearly, for a single layer graph/network, $\mathbf{e}_{\alpha}$ is a scalar
$1$ and the fourth-order tensor degenerates to the adjacency matrix of the
normal graphs. Similar to $A_{ij}$ in an adjacency matrix which indicates the
direction from the node $v_{j}$ to $v_{i}$, $A_{\alpha i\beta j}$ also
indicates the direction from the node $v_{\beta j}$ to the node $v_{\alpha i}$
in a network. Note that, vectors $\mathbf{e}_{i}$ and $\mathbf{e_{\alpha}}$
are not eigenvectors of the adjacency tensor. They are merely the vectors
characterizing features of the entities and layers, respectively. We shall
discuss the MLN-based spectral space in Section IV.

Given an adjacency tensor, we can define the Laplacian tensor of the
multilayer networks similar to that in a single-layer graph. Denoting the
degree (or multi-strength) of the entity’s $i$’s projected node $v_{\alpha i}$
in layer $\alpha$ as $d(v_{\alpha i})$ which is a summation over weights of
different natures (inter- and intra- layer edges), the degree tensor
$\mathbf{D}\in\mathbb{R}^{M\times N\times M\times N}$ is defined as a diagonal
tensor with entries $D_{\alpha i\alpha i}=d(v_{\alpha i})$ for $1\leq i\leq
N,1\leq\alpha\leq M$, whereas its other entries are zero. The Laplacian tensor
can be defined as follows.

######  Definition 3 (Laplacian Tensor).

A multilayer network $\mathcal{M}$, with $|\mathcal{L}|=M$ layers and
$|l_{i}|=N$ nodes in each layer $i$, can be represented by a fourth-order
Laplacian tensor $\mathbf{L}\in\mathbb{R}^{M\times N\times M\times N}$ defined
as $\mathbf{L=D-A}$, where $\mathbf{A}$ is the adjacency tensor and
$\mathbf{D}$ is the degree tensor.

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 MLN depending on datasets and user
objectives.

####  III-B2 Representation of General MLN

Representing a general multilayer network 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 [17] 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 MLN by the forth-order tensor in M-GSP.
Note that, although M-GSP also needs a reconstruction to represent a general
MLN, 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 spectral 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. (6), the
entities and layers can be virtual and not necessarily physical to capture the
underlying structures of the datasets. The information within the multilayer
networks, together with definitions of the underlying virtual entities and
layers, should only depend on the structure of the multilayer networks. We
will illustrate this further in Section VI-B.

###  III-C Flattening and Analysis

In this part, we introduce the flattening of the multilayer network, which
could simplify some operations in the tensor-based M-GSP. For a multilayer
network $\mathcal{M}=\\{\mathcal{V},\mathcal{L},\mathbf{F}\\}$ with $M$ layers
and $N$ nodes in each layer, its fourth-order representing tensor
$\mathbf{F}\in\mathbb{R}^{M\times N\times M\times 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 $\mathbf{F}$ can be flattened
into $\mathbf{F}_{FL}\in\mathbb{R}^{MN\times MN}$ with each element

${[F_{FL}]}_{{N(\alpha-1)+i,N(\beta-1)+j}}=F_{\alpha i\beta j}.$ (8)

  * •

Entity-wise Flattening: The representing tensor $\mathbf{F}$ can be flattened
into $\mathbf{F}_{FN}\in\mathbb{R}^{NM\times NM}$ with each element

${[F_{FN}]}_{{M(i-1)+\alpha,M(j-1)+\beta}}=F_{\alpha i\beta j}.$ (9)

These two flattening methods provide two ways to interpret the network
structure. In the first method, the flattened multilayer network 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 network 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 network with $3$ nodes in each layer are shown in
Fig. 3. From the examples, we can see that the diagonal blocks in
$\mathbb{R}^{N\times N}$ are the intralayer connections for each layer and
other blocks describe the interlayer connections through layer-wise
flattening; and the diagonal block in $\mathbb{R}^{M\times 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 MLN with a different indexing of vertices, they are still
helpful to analyze the MLN spectral space. For example, in [12], the
approximations of spectral radius are derived based on different structures of
these two flattened matrix.

Figure 3: Example of Multilayer Network Flattening.

###  III-D Signals and Shifting over the Multilayer Networks

Based on the tensor representation, we now define signals and signal shifting
over the multilayer networks. 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
network. Intuitively, if there are $K=MN$ nodes, there are $MN$ signal samples
in total. Similar to GSP, we use the representing (adjacency/Laplacian) tensor
$\mathbf{F}\in\mathbb{R}^{M\times N\times M\times N}$ as the basic MLN-filter.
Since the input signal and the output signal of the MLN-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.

######  Definition 4 (Signals over Multilayer Networks).

For a multilayer network
$\mathcal{M}=\\{\mathcal{V},\mathcal{L},\mathbf{F}\\}$, with $|\mathcal{L}|=M$
layers and $|l_{i}|=N$ nodes in each layer $i$, the definition of multilayer
network signals is a second-order tensor

$\mathbf{s}=(s_{\alpha i})\in\mathbb{R}^{M\times N},\quad 1\leq\alpha\leq
M,1\leq i\leq N,$ (10)

where the entry $s_{\alpha i}$ is the signal sample in the projected node of
entity $i$ in layer $\alpha$.

Note that, if the multilayer network degenerates to a single-layer graph with
$M=1$, the multilayer network signal becomes an $N$-length vector, which is
similar to that in the traditional GSP. Similar to the representing tensor,
the tensor signal $\mathbf{s}\in\mathbb{R}^{M\times N}$ can also be flattened
as a vector in $\mathbb{R}^{MN}$:

  * •

Layer-wise flattening: $\mathbf{s}_{L}\in\mathbb{R}^{MN}$ whose entries are
calculated as $[s_{L}]_{N(\alpha-1)+i}=s_{\alpha i}$.

  * •

Entity-wise flattening: $\mathbf{s}_{N}\in\mathbb{R}^{NM}$ whose entries are
calculated as $[s_{N}]_{M(i-1)+\alpha}=s_{\alpha i}$.

Given the definitions of multilayer network 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 network based on
the contraction (inner product) between the representing tensor and tensor
signals.

######  Definition 5 (Signal Shifting over Multilayer Networks).

Given the representing matrix $\mathbf{F}\in\mathbb{R}^{M\times N\times
M\times N}$ and the tensor signal $\mathbf{s}\in\mathbb{R}^{M\times N}$
defined over a multilayer network $\mathcal{M}$, the signal shifting is
defined as the contraction (inner product) between $\mathbf{F}$ and
$\mathbf{s}$ in one entity-related order and one layer-related order, i.e.,

$\mathbf{s}^{\prime}=\mathbf{F}\diamond\mathbf{s}\in\mathbb{R}^{M\times N},$
(11)

where $\diamond$ is the contraction between $\mathbf{F}$ and $\mathbf{s}$
defined in Eq. (4).

The elements in the shifted signal $\mathbf{s}^{\prime}$ are calculated as

$s^{\prime}_{\alpha i}=\sum_{\beta=1}^{M}\sum_{j=1}^{N}F_{\alpha i\beta
j}s_{\beta j}.$ (12)

From Eq. (12), there are two important factors to construct the shifted
signal: 1) The signal in the neighbors of the node $v_{\alpha i}$; and 2) The
intensity of interactions between the node $v_{\alpha i}$ and its neighbors.
Then, the signal shifting is related to the diffusion process over the
multilayer networks. More specifically, if $\mathbf{F}$ is the adjacency
tensor, signals shift in directions of edges. Meanwhile, if $\mathbf{F}$ is
the Laplacian tensor, Eq. (12) can be written as

$s^{\prime}_{\alpha i}=\sum_{\beta=1}^{M}\sum_{j=1}^{N}A_{\alpha i\beta
j}(s_{\alpha i}-s_{\beta j}),$ (13)

which is the weighted average of difference with neighbors.

##  IV Multilayer Network Spectral Space

In traditional GSP, graph spectral space is defined according to the
eigenspace of the representing matrix [3]. Similarly, we define the MLN
spectral 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
networks in this section for simplicity and clarity of presentation. For the
directed MLN, we provide alternative spectral definitions in Appendix C and
leave the frequency analysis in the future works. Meanwhile, all the proofs of
the properties in this part are listed in Appendix A.

###  IV-A Joint Spectral Analysis in M-GSP

For a multilayer network
$\mathcal{M}=\\{\mathcal{V},\mathcal{L},\mathbf{F}\\}$ with $M$ layers and $N$
nodes, the eigen-tensor $\mathbf{V}\in\mathbb{R}^{M\times N}$ of the
representing tensor $\mathbf{F}$ is defined in the tensor-based multilayer
network theory [11] as $\mathbf{F}\diamond\mathbf{V}=\lambda\mathbf{V}$. More
specifically, $\mathbf{F}\in\mathbb{R}^{M\times N\times M\times N}$ can be
decomposed as

$\displaystyle\mathbf{F}=\sum_{k=1}^{MN}\lambda_{k}\mathbf{V}_{k}\circ\mathbf{V}_{k}=\sum_{\alpha=1}^{M}\sum_{i=1}^{N}\lambda_{\alpha
i}\mathbf{V}_{\alpha i}\circ\mathbf{V}_{\alpha i},$ (14)

where $\lambda_{k}$ is the eigenvalues and
$\mathbf{V}_{k}\in\mathbb{R}^{M\times N}$ is the corresponding eigen-tensor.
Note that $\mathbf{V}_{\alpha i}$ just relabels the index of $\mathbf{V}_{k}$,
and there is no specific order for $\mathbf{V}_{\alpha 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 MLN Fourier space
as follows.

######  Definition 6 (Joint Multilayer Network Fourier Space).

For a multilayer network
$\mathcal{M}=\\{\mathcal{V},\mathcal{L},\mathbf{F}\\}$ with $M$ layers and $N$
nodes, the MLN Fourier space is defined as the space consisting of all
spectral tensors $\\{\mathbf{V}_{1},\cdots,\mathbf{V}_{MN}\\}$, which
characterizes the joint features of entities and layers.

Recall that in GSP, the GFT is defined based on the inner product of
$\mathbf{V}^{-1}$ and the signals $\mathbf{s}$ defined in Eq. (1). Similarly,
we can define the M-GFT based on the spectral tensors of the representing
tensor $\mathbf{F}$ to capture joint features of inter- and intra- layer
interactions as follows.

######  Definition 7 (Joint M-GFT).

Let $\mathbf{U}_{\mathcal{F}}=(\mathbf{V}_{\alpha i})\in\mathbb{R}^{M\times
N\times M\times N}$ consist of spectral tensors of the representing tensor
$\mathbf{F}$, where $[U_{\mathcal{F}}]_{\alpha i\beta j}=[V_{\alpha i}]_{\beta
j}$.

The joint M-GFT can be defined as the contraction between
$\mathbf{U}_{\mathcal{F}}$ and the tensor signal
$\mathbf{s}\in\mathbb{R}^{M\times N}$, i.e.,

$\hat{\mathbf{s}}=\mathbf{U}_{\mathcal{F}}\diamond\mathbf{s}.$ (15)

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. (8) and Eq. (9) lead to the same
eigenvalues.

This property shows that the flattened tensors are the reshape of the original
representing tensor, and could capture some of the spectral properties as
follows.

######  Property 2.

Given the eigenpair $(\lambda_{FL},\mathbf{x})$ of the layer-wise flattened
tensor, the eigenpair $(\lambda,\mathbf{V})$ of the original representing
tensor can be calculated as $\lambda=\lambda_{FL}$, and $V_{\alpha
i}=x_{N(\alpha-1)+i}$. Similarly, given the eigenpair
$(\lambda_{FN},\mathbf{y})$ of the entity-wise flattened tensor, the eigenpair
$(\lambda,\mathbf{V})$ of the original representing tensor can be calculated
as $\lambda=\lambda_{FN}$, and $V_{\alpha i}=y_{M(i-1)+\alpha}$.

The Property 2 shows that we can calculate the eigen-tensor from the flattened
tensor to simplify the decomposition operations. Moreover, the joint M-GFT is
the bijection of GFT in the flattened MLN, with vertices indexed by both the
layers and the entities. However, such joint M-GFT 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.

###  IV-B Order-wise Spectral Analysis in M-GSP

In an undirected multilayer network, the representing tensor
(adjacency/Laplacian) $\mathbf{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
network structure under orthogonal CP-decomposition [26] as follows:

$\displaystyle\mathbf{F}$
$\displaystyle\approx\sum_{\alpha=1}^{M}\sum_{i=1}^{N}\lambda_{\alpha
i}\cdot\mathbf{f}_{\alpha}\circ\mathbf{e}_{i}\circ\mathbf{f}_{\alpha}\circ\mathbf{e}_{i}$
(16) $\displaystyle=\sum_{\alpha=1}^{M}\sum_{i=1}^{N}\lambda_{\alpha
i}\tilde{\mathbf{V}}_{\alpha i}\circ\tilde{\mathbf{V}}_{\alpha i},$ (17)

where $\mathbf{f}_{\alpha}\in\mathbb{R}^{M}$ are orthonormal,
$\mathbf{e}_{i}\in\mathbb{R}^{N}$ are orthonormal and
$\tilde{\mathbf{V}}_{\alpha
i}=\mathbf{f}_{\alpha}\circ\mathbf{e}_{i}\in\mathbb{R}^{M\times N}$.

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 Appendix B. In Eq. (16), $\mathbf{f}_{\alpha}$ and
$\mathbf{e}_{i}$ capture the features of layers and entities, respectively,
which can be interpreted as the subspaces of the MLN. More discussions about
the frequency interpretation of order-wise M-GSP spectrum and connections to
MWGSP spectrum are presented in Section VI-C.

Note that, if there is only one layer in the multilayer network, Eq. (16)
reduces to the eigendecomposition of a normal single-layer graph, i.e.,
$\mathbf{F}=\sum_{i=1}^{N}\lambda_{i}\mathbf{e}_{i}\circ\mathbf{e}_{i}$

With the decomposed representing tensor in Eq. (16), the order-wise MLN
spectrum is defined as follows.

######  Definition 8 (Order-wise MLN Spectral Pair).

For a multilayer network
$\mathcal{M}=\\{\mathcal{V},\mathcal{L},\mathbf{F}\\}$ with $M$ layers and $N$
nodes, the order-wise MLN spectral pairs are defined by $\\{\lambda_{\alpha
i},\mathbf{f}_{\alpha},\mathbf{e}_{i}\\}$, where
$\\{\mathbf{f}_{1},\cdots,\mathbf{f}_{M}\\}$ and
$\\{\mathbf{e}_{1},\cdots,\mathbf{e}_{N}\\}$ characterize features of layers
and entities, respectively.

With the definition of order-wise MLN spectral pair, we now explore their
properties. Considering $\tilde{\mathbf{V}}_{\alpha
i}=\mathbf{f}_{\alpha}\circ\mathbf{e}_{i}$, we have the following property,
which indicates the availability of a joint MLN analysis based on order-wise
spectrum.

######  Property 3.

The factor tensor $\tilde{\mathbf{V}}_{\alpha i}$ of the representing tensor
$\mathbf{F}$ is the approximated eigen-tensor of $\mathbf{F}$.

By constructing a fourth-order tensor
$\tilde{\mathbf{U}}_{\mathcal{F}}\in\mathbb{R}^{M\times N\times M\times N}$
with $\tilde{\mathbf{V}}_{\alpha i}$ as its elements, i.e.,
$[\tilde{U}_{\mathcal{F}}]_{\alpha i\beta j}=[\tilde{V}_{\alpha i}]_{\beta
j}$, we can have the following property.

######  Property 4.

Let
$\mathbf{W}=\tilde{\mathbf{U}}_{\mathcal{F}}\otimes\tilde{\mathbf{U}}_{\mathcal{F}}$,
where $\otimes$ is the contraction in the third and forth order with
$W_{\alpha i\beta j}=\sum_{p,\theta}[\tilde{U}_{\mathcal{F}}]_{\beta j\theta
p}\times[\tilde{U}_{\mathcal{F}}]_{\alpha i\theta p}$. Then, $\mathbf{W}$ is
super-diagonal with super-diagonal elements all equal to one.

This property generalizes the orthogonality of the spectral tensor into a
similar definition of matrix.

We now introduce the order-wise MLN spectral transform. Similar to Eq. (15),
the joint transform can be defined as

$\hat{\mathbf{s}}=\tilde{\mathbf{U}}_{\mathcal{F}}\diamond\mathbf{s}.$ (18)

Note that each element of $\hat{\mathbf{s}}$ in Eq. (18) can be calculated as

$\displaystyle\hat{s}_{\alpha i}$
$\displaystyle=\sum_{\beta,j}[\tilde{U}_{\mathcal{F}}]_{\alpha i\beta
j}s_{\beta j}=\sum_{\beta,j}[\tilde{V}_{\alpha i}]_{\beta j}s_{\beta j}$ (19)
$\displaystyle=\sum_{\beta,j}[f_{\alpha}]_{\beta}\cdot[e_{i}]_{j}\cdot
s_{\beta j}.$ (20)

Let $\mathbf{E}_{f}=[\mathbf{f}_{1}\cdots\mathbf{f}_{M}]\in\mathbb{R}^{M\times
M}$ and
$\mathbf{E}_{e}=[\mathbf{e}_{1}\cdots\mathbf{e}_{N}]\in\mathbb{R}^{N\times
N}$. We then have
$\hat{\mathbf{s}}^{\prime}=\mathbf{E}_{f}^{\mathrm{T}}\mathbf{s}\mathbf{E}_{e}$,
with each element
$\tilde{s}_{\alpha}^{\prime}=\sum_{j,\beta}[f_{\alpha}]_{\beta}\cdot[e_{i}]_{j}\cdot
s_{\beta j}$. Clearly, the joint M-GFT can be obtained as
$\hat{\mathbf{s}}=\hat{\mathbf{s}}^{\prime}=\mathbf{E}_{f}^{\mathrm{T}}\mathbf{s}\mathbf{E}_{e}$.
Then, we have the following definition of M-GFT based on order-wise spectrum.

######  Definition 9 (Order-wise M-GFT).

Given the spectral vectors
$\mathbf{E}_{f}=[\mathbf{f}_{1}\cdots\mathbf{f}_{M}]\in\mathbb{R}^{M\times M}$
and $\mathbf{E}_{e}=[\mathbf{e}_{1}\cdots\mathbf{e}_{N}]\in\mathbb{R}^{N\times
N}$, the layer-wise M-GFT can be defined as

$\hat{\mathbf{s}}_{L}=\mathbf{E}_{f}^{\mathrm{T}}\mathbf{s}\in\mathbb{R}^{M\times
N},$ (21)

and the entity-wise M-GFT can be defined as

$\hat{\mathbf{s}}_{N}=\mathbf{s}\mathbf{E}_{e}\in\mathbb{R}^{M\times N}.$ (22)

The joint M-GFT based on order-wise spectrum is defined as

$\hat{\mathbf{s}}=\mathbf{E}_{f}^{\mathrm{T}}\mathbf{s}\mathbf{E}_{e}\in\mathbb{R}^{M\times
N}.$ (23)

If there is only one layer in the multilayer network, the M-GFT calculated
with $\mathbf{s}^{T}\in\mathbb{R}^{N}$ as
$(\hat{\mathbf{s}}_{N})^{\mathrm{T}}=(\mathbf{s}\mathbf{E}_{e})^{\mathrm{T}}\in\mathbb{R}^{N}$,
which has the same form as the traditional GFT in Eq. (1).

In addition, since $\mathbf{f}_{\alpha}$ and $\mathbf{e}_{i}$ are orthonormal
basis of undirected MLN, the inverse M-GFT can be calculated as

$\mathbf{s}^{\prime}=\mathbf{E}_{f}\hat{\mathbf{s}}\mathbf{E}_{e}^{\mathrm{T}}.$
(24)

Different from joint MLN Fourier space in Section IV-A, the order-wise MLN
spectrum provides an individual analysis on layers and entities separately,
and a reliable approximated analysis on the underlying MLN structures jointly.

###  IV-C MLN 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 networks.

Given the multilayer network
$\mathcal{M}=\\{\mathcal{V},\mathcal{L},\mathbf{F}\\}$ with $M$ layers and $N$
nodes in each layer, its representing tensor $\mathbf{F}\in\mathbb{R}^{M\times
N\times M\times N}$ can be decomposed via HOSVD as

$\mathbf{F}=\mathbf{S}\times_{1}\mathbf{U}^{(1)}\times_{2}\mathbf{U}^{(2)}\times_{3}\mathbf{U}^{(3)}\times_{4}\mathbf{U}^{(4)},$
(25)

where
$\mathbf{U}^{(n)}=[\mathbf{U}^{(n)}_{1}\quad\mathbf{U}^{(n)}_{2}\quad\cdots\quad\mathbf{U}^{(n)}_{I_{n}}]$
is a unitary $(I_{n}\times I_{n})$ matrix, with $I_{1}=I_{3}=M$ and
$I_{2}=I_{4}=N$. $\mathbf{S}$ is a complex $(I_{1}\times I_{2}\times
I_{3}\times I_{4})$-tensor of which the subtensor $\mathbf{S}_{i_{n}}$
obtained by fixing $n$th index to $\alpha$ have

  * •

$<\mathbf{S}_{i_{n}=\alpha},\mathbf{S}_{i_{n}=\beta}>=0$ where
$\alpha\neq\beta$.

  * •

$||\mathbf{S}_{i_{n}=1}||\geq||\mathbf{S}_{i_{n}=2}||\geq\cdots\geq||\mathbf{S}_{i_{n}=I_{n}}||\geq
0$.

The Frobenius-norms $\sigma_{i}^{(n)}=||S_{i_{n}=i}||$ is the $n$-mode
singular value, and $\mathbf{U}^{(i)}$ are the corresponding $n$-mode singular
vectors. For an undirected multilayer network, the representing tensor is
symmetric for every 2-D combination. Thus, there are two modes of singular
spectrum, i.e., $(\gamma_{\alpha},\mathbf{f}_{\alpha})$ for mode $1,3$, and
$(\sigma_{i},\mathbf{e}_{i})$ for mode $2,4$. More specifically,
$\mathbf{U}^{(1)}=\mathbf{U}^{(3)}=(\mathbf{f}_{\alpha})$ and
$\mathbf{U}^{(2)}=\mathbf{U}^{(4)}=(\mathbf{e}_{i})$. Since the joint singular
tensor captures the consistent information of entities and layers, it can be
calculated as

$(\lambda_{\alpha i},\hat{\mathbf{V}}_{\alpha
i})=(\gamma_{\alpha}\cdot\sigma_{i},\mathbf{f}_{\alpha}\cdot\mathbf{e}_{i}).$
(26)

Note that the diagonal entries of $\mathbf{S}$ are not the eigenvalues or
frequency coefficients of the representing tensor in general. The multilayer
network singular space is defined as follows.

######  Definition 10 (Multilayer Network Singular Space).

For a multilayer network
$\mathcal{M}=\\{\mathcal{V},\mathcal{L},\mathbf{F}\\}$ with $M$ layers and $N$
nodes, the MLN singular space is defined as the space consisting of all
singular tensors $\\{\hat{\mathbf{V}}_{1}\cdots\hat{\mathbf{V}}_{MN}\\}$
obtained from Eq. (26). The singular vectors
$\\{\mathbf{f}_{1},\cdots,\mathbf{f}_{M}\\}$ and
$\\{\mathbf{e}_{1},\cdots,\mathbf{e}_{N}\\}$ in Eq. (25) characterize layers
and entities, respectively.

Similar to order-wise spectral analysis in Section IV-B, we can define the MLN
singular tensor transform (M-GST) based on the singular tensors as follows.

######  Definition 11 (M-GST).

Suppose that
$\mathbf{U}_{s}=(\mathbf{f}_{\alpha}\circ\mathbf{e}_{i})\in\mathbb{R}^{M\times
N\times M\times N}$ consists of the singular vectors of the representing
tensor $\mathbf{F}$ in Eq. (25), where $[U_{s}]_{\alpha i\beta
j}=[f_{\alpha}]_{\beta}\cdot[e_{i}]_{j}$. The M-GST can be defined as the
contraction between $\mathbf{U}_{s}$ and the tensor signal
$\mathbf{s}\in\mathbb{R}^{M\times N}$, i.e.,

$\check{\mathbf{s}}=\mathbf{U}_{s}\diamond\mathbf{s}.$ (27)

If the singular vectors are included in
$\mathbf{W}_{f}=[\mathbf{f}_{1}\cdots\mathbf{f}_{M}]\in\mathbb{R}^{M\times M}$
and $\mathbf{W}_{e}=[\mathbf{e}_{1}\cdots\mathbf{e}_{N}]\in\mathbb{R}^{N\times
N}$, the layer-wise M-GST can be defined as

$\check{\mathbf{s}}_{L}=\mathbf{W}_{f}^{\mathrm{T}}\mathbf{s}\in\mathbb{R}^{M}\times\mathbb{R}^{N},$
(28)

and the entity-wise M-GST can be defined as

$\check{\mathbf{s}}_{N}=\mathbf{s}\mathbf{W}_{e}\in\mathbb{R}^{M}\times\mathbb{R}^{N}.$
(29)

Inverse M-GST can be defined similarly as in Eq. (24) with unitary
$\mathbf{W}_{e}$ and $\mathbf{W}_{f}$.

Compared to the eigen-tensors in Eq. (14), 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, together with
additional discussions in Appendix B. The intuition of applying HOSVD in MLN
analysis and its correlations to GSP are also provided in Section VI-A3.

###  IV-D Spectrum Ranking in the Multilayer Network

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
$\lambda_{\alpha i}$ based on the spectral tensors to rank the MLN
frequencies. Let $|\lambda|_{max}$ be the joint singular/eigen- value with the
largest magnitude. The M-GSP total variation is defined as follows:

$\displaystyle TV(\mathbf{V}_{\alpha i})$ $\displaystyle=||\mathbf{V}_{\alpha
i}-\frac{1}{|\lambda|_{max}}\mathbf{F}\diamond\mathbf{V}_{\alpha i}||_{1}$
(30)
$\displaystyle=|1-\frac{\lambda}{|\lambda|_{max}}\mathbf{|}\cdot||{V}_{\alpha
i}||_{1},$ (31)

where $||\cdot||_{1}$ is the $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 $\lambda_{\alpha i}$
is said to be a higher frequency if its total variation $TV(\mathbf{V}_{\alpha
i})$ is larger, and its corresponding spectral tensor $\mathbf{V}_{\alpha i}$
is a higher frequency spectrum. We shall provide more details on
interpretation of MLN frequency in Section VI-A.

##  V Filter Design

In this section, we introduce an MLN filter design together with its
properties based on signal shifting.

###  V-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.,
$\mathbf{s}^{\prime}=f_{1}(\mathbf{s})=\mathbf{F}\diamond\mathbf{s}$.
Similarly, a second order filter can be defined as additional filtering on
first-order filtered signal, i.e.,

$\displaystyle\mathbf{s}^{\prime\prime}=f_{2}(\mathbf{s})=\mathbf{F}\diamond(\mathbf{F}\diamond\mathbf{s}),$
(32)

whose entries $s_{\alpha i}^{\prime\prime}$ are calculated as

$\displaystyle s_{\alpha i}^{\prime\prime}$
$\displaystyle=\sum_{\beta=1}^{M}\sum_{j=1}^{N}F_{\alpha i\beta j}s_{\beta
j}^{\prime}=\sum_{\beta,j}F_{\alpha i\beta j}\sum_{\epsilon,p}F_{\beta
j\epsilon p}s_{\epsilon p}$ (33) $\displaystyle=\sum_{\epsilon,p}s_{\epsilon
p}\sum_{\beta,j}F_{\alpha i\beta j}F_{\beta j\epsilon
p}=(\mathbf{F}\odot\mathbf{F})\diamond\mathbf{s},$ (34)

where $\odot$ is the contraction defined in Eq. (5).

Let $\mathbf{F}^{[2]}=\mathbf{F}\odot\mathbf{F}$. From Eq. (14), we have:

$\displaystyle F^{[2]}_{\alpha i\beta j}$
$\displaystyle=\sum_{\theta,p}F_{\alpha i\theta p}F_{\theta p\beta j}$ (35)
$\displaystyle=\sum_{\theta,p}(\sum_{k}\lambda_{k}[{V}_{k}]_{\alpha
i}[V_{k}]_{\theta p})(\sum_{t}\lambda_{t}[V_{t}]_{\beta j}[V_{t}]_{\theta p})$
$\displaystyle=\sum_{k,t}\lambda_{k}\lambda_{t}[{V}_{k}]_{\alpha
i}[{V}_{t}]_{\beta j}(\sum_{\theta,p}[{V}_{t}]_{\theta p}[{V}_{k}]_{\theta
p})$ $\displaystyle=\sum_{k}\lambda_{k}^{2}[{V}_{k}]_{\alpha
i}[{V}_{k}]_{\beta j}.$ (36)

Similarly, for $\tau$th-order term $\mathbf{F}^{[\tau]}$, its entry
$F^{[\tau]}_{\alpha i\beta j}$ can be calculated as $F^{[\tau]}_{\alpha i\beta
j}=\sum_{k}\lambda_{k}^{\tau}[\mathbf{V}_{k}]_{\alpha
i}[\mathbf{V}_{k}]_{\beta j}$.

Now we have the following property.

######  Property 5.

The $\tau$-th order basic shifting filter $f_{\tau}(\mathbf{s})$ can be
calculated as

$\displaystyle f_{\tau}(\mathbf{s})$
$\displaystyle=\mathbf{F}^{[\tau]}\diamond\mathbf{s}=(\sum_{k=1}^{MN}\lambda_{k}^{\tau}\mathbf{V}_{k}\circ\mathbf{V}_{k})\diamond\mathbf{s}.$
(37)

This property is the M-GSP counterpart to the traditional linear system
interpretation that complex exponential signals are eigenfunction of linear
systems [3], and provides a quicker implementation of higher-order shifting.
With the $k$-order polynomial term, the adaptive polynomial filter is defined
as

$h(\mathbf{s})=\sum_{k}\alpha_{k}\mathbf{F}^{[k]}\diamond\mathbf{s},$ (38)

where $\\{\alpha_{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 MLN polynomial filtering in Section VII.

###  V-B Spectral Filter Design

Filtering in the graph spectral space is useful in GSP frequency analysis. For
example, ordering the Laplacian graph spectrum
$\mathbf{V}_{\mathcal{G}}=[\mathbf{e}_{1},\cdots,\mathbf{e}_{N}]\in\mathbb{R}^{N\times
N}$ in a descent order by the graph total variation [3], i.e., high frequency
to low frequency, the GFT of $\mathbf{s}\in\mathbb{R}^{N}$ is calculated as
$\hat{\mathbf{s}}=\mathbf{V}_{\mathcal{G}}^{\mathrm{T}}\mathbf{s}$. By
removing $k$ elements in the low frequency part, i.e.,
$\hat{s}^{\prime}=[\hat{s}_{1},\cdots,\hat{s}_{N-k},0,\cdots,0]$, a high-pass
filter can be designed as

$\displaystyle\mathbf{s}^{\prime}$
$\displaystyle=\mathbf{V}_{\mathcal{G}}\hat{s}^{\prime}=\mathbf{V}_{\mathcal{G}}\Sigma_{k}\mathbf{V}_{\mathcal{G}}^{\mathrm{T}}\mathbf{s}$
(39)

where $\Sigma_{k}=diag([\sigma_{1},\cdots,\sigma_{N}])$ is a diagonal matrix
with $\sigma_{i}=0$ for $i=1,\cdots,N-k$; otherwise, $\sigma_{i}=0$.

Similarly, in M-GSP, a spectral filter is designed by filtering in the
spectral space together with the inverse M-GFT. With Eq. (23) and Eq. (27),
spectral filtering of $\mathbf{s}$ is defined as

$\displaystyle\mathbf{s}^{\prime}=$
$\displaystyle\mathbf{E}_{f}\begin{bmatrix}g(\gamma_{1})&\cdots&0\\\
\vdots&\ddots&\vdots\\\
0&\cdots&g(\gamma_{N})\end{bmatrix}\mathbf{E}_{f}^{\mathrm{T}}\mathbf{s}\mathbf{E}_{e}\begin{bmatrix}f(\sigma_{1})&\cdots&0\\\
\vdots&\ddots&\vdots\\\
0&\cdots&f(\sigma_{N})\end{bmatrix}\mathbf{E}_{e}^{\mathrm{T}}$ (40)

where functions $g(\cdot)$ and $f(\cdot)$ are designed by the specific tasks.
For example, if one wants to design a layer-wise filter capturing the
smoothness of signals in the MLN singular space, the function $g(\cdot)$ can
be designed as $\mathbf{g}=[1,\cdots,1,0,\cdots,0]$ by ordering the layer-wise
singular vectors in the descent order of singular values. More discussions and
examples are presented in Section VI and Section VII.

In addition to polynomial and spectral filters, filters designed through
optimization based on geometric information play an important role in semi-
supervised signal processing works. Interested readers can refer to Appendix E
for a short discussion.

##  VI Discussion and Interpretative Insights

###  VI-A Interpretation of M-GSP Frequency

####  VI-A1 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] or Appendix D 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.,
$\mathbf{s}=[s_{1},\cdots,s_{N}]\in\mathbb{R}^{N}$, the signal shifting can be
interpreted as the shift filtering corresponding to a cyclic graph shown in
Fig. 4. Suppose that its adjacency matrix is written as

$\displaystyle{\mathbf{C}_{N}=\left[\begin{smallmatrix}0&0&\cdots&0&1\\\
1&0&\cdots&0&0\\\ \vdots&\ddots&\ddots&\ddots&\vdots\\\ 0&0&\ddots&0&0\\\
0&0&\cdots&1&0\end{smallmatrix}\right]}$ (41)

Then, the shifted signal over the cyclic graph is calculated as
$\mathbf{s}^{\prime}=\mathbf{C}_{N}\mathbf{s}=[s_{N}\quad
s_{1}\quad\cdots\quad s_{N-1}]^{\mathrm{T}}$, which shifts the signal at each
node to its next node. More specifically, $\mathbf{C}_{N}$ can be decomposed
as $\mathbf{C}_{N}=\mathbf{V\Lambda}\mathbf{V}^{-1}$ where the eigenvalues
$\lambda_{n}=e^{-j\frac{2\pi n}{N}}$ and
$\mathbf{V}^{-1}=\frac{1}{\sqrt{N}}[\lambda^{kn}_{N}]$ is the discrete Fourier
matrix. Inspired by the DSP, the eigenvectors in $\mathbf{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
(Laplacian/adjacency) matrix $\mathbf{F}_{M}$ of an arbitrary graph, the graph
Fourier space consists of the eigenvectors of $\mathbf{F}_{M}$ and the graph
frequencies are related to the eigenvalues. More specifically, the graph
Fourier space can be interpreted as the manifold or spectral 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.,
$TV(\mathbf{e}_{i})=||\mathbf{e}_{i}-\frac{1}{|\lambda|_{max}}\mathbf{e}_{i}||_{1}$,
to rank the graph frequencies, where $\mathbf{e}_{i}$ is the spectral
component related to the eigenvalue $\lambda_{i}$ in $\mathbf{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.

Figure 4: Example of Cyclic Graph.

####  VI-A2 Interpretation of MLN Frequency

Now, return to M-GSP. Given spectral tensors
$\mathbf{V}_{k}\in\mathbb{R}^{M\times N}$ of a multilayer network, a signal
$\mathbf{s}\in\mathbb{R}^{M\times N}$ can be written in a weighted sum of the
spectrum, i.e., $\mathbf{s}=\sum_{k}a_{k}\mathbf{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. (30).
If the signal component changes faster over the multilayer network, the
corresponding total variation is larger. Since we relate higher frequency
component with a larger total variation, the MLN frequency indicates how fast
the signal propagates over the multilayer network under the representing
tensor. If a signal $\mathbf{s}$ contains more high frequency components, it
changes faster under the representing matrix.

Figure 5: Example of MLN Frequencies. TABLE I: Total Difference of all Nodes Eigenvalue | 0.5726 | 2.7500 | 5.7259
---|---|---|---
Total Differences | 3.7439 | 2.6846 | 0.8293

Here, we use an example to illustrate it further. We randomly generate a
multiplex network with six layers and five nodes in each layer. Each node has
a probability of $30\%$ to connect other nodes for intralayer connections and
its counterparts in other layers for interlayer connections. We use the
Laplacian tensor as representing tensor and $l_{1}$-norm of the flattened
signal to calculate the total variation. Then we have $\lambda_{k}\geq 0$ and
$TV(\mathbf{V}_{k})=|1-\frac{\lambda}{\lambda_{max}}\mathbf{|}$ for $1\leq
k\leq MN$. Clearly, the total variation is larger if the eigenvalue is
smaller, i.e., smaller eigenvalues correspond to higher frequencies, similar
to GSP [3]. We next evaluate the differences between filtered signal and
original signal by treating different eigen-tensor as signals. From results of
Fig. 5, shifted signal sample of each node changes more than its original
samples. Also from the results given n Table I, we can see that a higher
frequency signal component exhibits a larger total difference between itself
and its shifts, indicating larger oscillations and faster propagation over the
MLN under its representing tensor.

####  VI-A3 Interpretation of MLN Singular Tensors

As discussed in Section VI-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 MLN-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 MLN singular space (M-GST) for the analysis of underlying
structures can be a useful alternative for M-GFT.

###  VI-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 networks
with different numbers of nodes, we may insert some isolated artificial nodes
to augment the multilayer network. 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
network. 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. 6, 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.

Figure 6: Example of Different Entities.

###  VI-C Distinctions from Existing GSP Works

####  VI-C1 Graph Signal Processing

Generally, M-GSP extends traditional GSP into multilayer networks. Although
one can stack all MLN 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 MLN, 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 intra- layer connections
together, while the M-GSP provide a more flexible choice on joint and order-
wise analysis. In Section IV-A, the joint M-GSP analysis introduced can be
viewed as the bijection of GFT in the flattened MLN, 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 MLN-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.

####  VI-C2 Joint Time-Vertex Fourier Analysis

In [17], a joint time-vertex Fourier transform (JFT) is defined by
implementing GFT and DFT consecutively. As discussed in Section VI-A1, the
time sequences can be interpreted under a cyclic graph, and thus reside on a
MLN structure. However, JFT assumes that all the layers have the same intra-
layer connections, which limits the generalization of the MLN analysis.
Differently, the tensor-based representation allows heterogeneous structures
for the intra-layer connections, which makes M-GSP more general.

####  VI-C3 Multi-way Graph Signal Processing

In [18], MWGSP has been proposed to process high-dimensional signals. Given
$K$th-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 network
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 MLN. Second, MWGSP assumes all the layers residing on a homogeneous
factor graph and restricts the types of manageable MLN 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 MLN. 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.

(a) Original Image.

(b) Superpixel.

(c) K-Means.

(d) GSP.

(e) MLN-FZ.

(f) MLN-SVD.

Figure 7: Example of BCCD Datasets and Segmented Images: (a) the original
image; (b) the boundaries of each superpixel; (c)-(f) segmented images under
different methods (WBCs are marked yellow, RBCs are marked green, and Platelet
(P) is marked blue).

##  VII Application Examples

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

###  VII-A Image Segmentation

In this part, we introduce an MLN spectral clustering for unsupervised RGB
image segmentation.

To model an RGB image using MLN, 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 multilayer network, as shown in Fig. 7(b). Here, we define
the feature of a superpixel according to its RGB pixel values. For interlayer
connections, each node connects with it counterparts in other layers. For
intralayer connections in layer $\ell$, we calculate the Gaussian-based
distance between two superpixels according to
$W_{ij}=\exp\left(-\frac{|\mathbf{s}_{i}(\ell)-\mathbf{s}_{j}(\ell)|^{2}}{\delta_{\ell}^{2}}\right)$
if $\quad|\mathbf{s}_{i}(\ell)-\mathbf{s}_{j}(\ell)|^{2}\leq t_{\ell}$;
otherwise, $W_{ij}=0$, where $\mathbf{s}_{i}(\ell)$ is the superpixel value in
layer $\ell$, $\delta_{\ell}$ is an adjustable parameter and $t_{\ell}$ 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 network with $M=3$ and $N$ nodes.

We now consider MLN-based spectral clustering. For image segmentation, we
focus on the properties of entities (i.e., superpixels), and implement
spectral clustering over entity-wise 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. (16), and singular vector in Eq. (25). 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.

TABLE II: Results of Image Segmentation in Image BSD300: (F) for all, and (C) for coarse | N=100 (F) | N=300 (F) | N=100 (C) | N=300 (C) | N=900 (C)
---|---|---|---|---|---
GSP | 0.1237 | 0.1149 | 0.3225 | 0.3087 | 0.3067
K-MEANS | 0.1293 | 0.1252 | 0.3044 | 0.3105 | 0.3124
MLN-SVD | 0.1326 | 0.1366 | 0.3344 | 0.3394 | 0.3335
MLN-CP | 0.1321 | 0.1293 | 0.3195 | 0.3256 | 0.3243

1: Input: RGB Image $\mathbf{I}\in\mathbb{R}^{P\times Q\times 3}$;

2: Build $N$ superpixels for the image $\mathbf{I}$ and calculate the value of
superpixel based on the mean of all pixels inside that superpixel, i.e.,
$\mathbf{s}\in\mathbb{R}^{1000\times 3}$;

3: Construct a multilayer network $\mathbf{A}\in\mathbb{R}^{M\times N\times
M\times N}$;

4: Find entity-wise spectrum
$\mathbf{E}=[\mathbf{e}_{1},\cdots,\mathbf{e}_{N}]\in\mathbb{R}^{N\times N}$;

5: Select the first $K$ important leading spectrum based on the eigenvalues
(singular values) of $\mathbf{E}$ as $\mathbf{C}\in\mathbb{R}^{N\times K}$;

6: Cluster each row of $\mathbf{C}$, and assign the $i$th superpixel into
$j$th cluster if the $i$th row of $\mathbf{C}$ is clustered into $j$th group;

7: Assign all pixels inside one superpixel to the cluster of that superpixel;

8: Output: The segmented image.

Algorithm 1 MLN-based Unsupervised Image Segmentation

To test the proposed Algorithm 1, we first compare its results with those from
GSP-based method and traditional $k$-means clustering by using a public BCCD
blood cell dataset shown in Fig. 7(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
$\sum_{\ell=1}^{3}|\mathbf{s}_{i}(\ell)-\mathbf{s}_{j}(\ell)|^{2}$ to form
edge connections in a single graph. There is only a single $\delta_{\ell}$ and
$t_{\ell}$ in the Gaussian model. For M-GSP, we use the MLN singular vectors
(MLN-SVD), and tensor factorization (MLN-FZ) for spectral clustering,
separately. Their respective results are compared in Fig 7(c)-7(f). WBCs are
marked yellow, and RBCs are marked green. Platelet (P) is marked blue. From
the illustrative results, MLN methods exhibit better robustness and are better
in detecting regions under noise. Comparing results from different MLN-based
methods, we find MLN-FZ to be less stable than HOSVD, partly due to
approximation algorithms used for tensor factorization. Overall, MLN-based
methods shows reliable and robust performance over GSP-based method and
$k$-means. Additional visualized results of BCCD images can be found in
Appendix F.

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
datasets [39]. 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 $J(A,B)=\frac{|A\cap
B|}{|A\cup B|}$. 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 and GSP-based spectral clustering.
The best performance is marked in bold. From the results of Table II, we can
see that larger number clusters of 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, MLN-based methods still demonstrate
better performance. Moreover, even when we use the same total number of nodes
in a single layer graph and multilayer network for another fairness comparison
in terms of complexity, i.e., $N=300$ for graph and $N=100$ for MLN, MLN-based
methods still perform better than graph-based methods in this example
application.

###  VII-B Semi-Supervised Classification

Semi-supervised classification is an important practical application for
graph-based signal processing. In this application, we apply MLN polynomial
filters for semi-supervised classification. Traditional GSP defines adaptive
filter as $f(\mathbf{s})=\sum_{i}a_{i}\mathbf{W}^{i}\mathbf{s}$, where
$\mathbf{W}$ is an adjacency matrix based on pairwise distance or a
representing matrix constructed from the adjacency matrix. Here, signals are
defined as labels or confidence values of nodes, i.e.,
$\mathbf{s}=[\mathbf{s}_{L}^{\mathrm{T}}\>\mathbf{0}_{UL}^{\mathrm{T}}]^{\mathrm{T}}$
by setting unlabeled signals to zero. To estimate parameters $a_{i}$ of
$f(\cdot)$, Optimization can be formulated to minimize, e.g., the mean square
error (MSE) from ground truth label $\mathbf{y}_{L}$

$\min_{\mathbf{a}}\quad||M(f(\mathbf{s}))_{L}-\mathbf{y}_{L}||_{2}^{2},$ (42)

where $M(\cdot)$ is a mapping of filtered signals to their discrete labels.
For example, in a $\\{\pm 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.

Figure 8: Results of Classification.

Similarly, in an MLN, we can also apply polynomial filters for label
estimation. Given an arbitrary dataset
$\mathbf{X}=[\mathbf{x}_{1},\cdots,\mathbf{x}_{N}]\in\mathbb{R}^{K\times N}$
with $N$ signal points and $K$ features for each node, we can construct a MLN
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
$\mathbf{A}\in\mathbb{R}^{M\times N\times M\times N}$. A signal is defined by

$\mathbf{s}=\begin{bmatrix}\mathbf{s}_{L}&\cdots&\mathbf{s}_{L}\\\
\mathbf{0}_{UL}&\cdots&\mathbf{0}_{UL}\end{bmatrix}^{\mathrm{T}}\in\mathbb{R}^{M\times
N},$ (43)

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. (43) for notational convenience. We now apply polynomial filters on
signals, i.e.,
$\mathbf{s}_{1}=h_{1}(\mathbf{s})=\sum_{i}a_{i}\mathbf{A}^{[i]}\diamond\mathbf{s}$,
and $\mathbf{s}_{2}=h_{2}(\mathbf{s})=\mathbf{A}^{[i]}\diamond\mathbf{s}$, For
a filtered signal $\mathbf{s}_{X}\in\mathbb{R}^{M\times N}$ ($X=1,2$), we
define a function to map 2-D signals into 1-D by calculated the column-wise
mean of $\mathbf{s}_{X}$, i.e., $\mathbf{\bar{s}}_{X}={\rm
mean}_{col}(\mathbf{s}_{X})\in\mathbb{R}^{1\times N}$.

Next, we can define a function $M(\cdot)$ on $\mathbf{\bar{s}}_{X}$ and
consider certain objective functions in filter design. To validate the
efficacy of polynomial filtering in the MLN framework, we test $h_{1}(\cdot)$
and $h_{2}(\cdot)$ 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 MLN with $N=297$ nodes in each of the
$M=13$ layers. More specifically, we directly use the labels as $\mathbf{s}$.
For $h_{1}(\cdot)$ (AF), we set $a_{i}\neq 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}(\cdot)$ (APF), we estimate a classification threshold via the mean of
$\bar{\mathbf{s}}_{X}$ 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
$\bar{\mathbf{s}}_{X}$ in M-GSP and use
$\mathbf{s}^{\prime}=f([\mathbf{s}_{L}^{T}\;\mathbf{0}_{\rm UL}^{T}]^{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. 8, GSP-based and M-GSP
based methods exhibit better performance than traditional learning algorithms,
particularly when the fraction of training samples is small. In general, M-GSP
based methods demonstrate superior performance among all methods owing to its
strength to extract ‘multilayer’ features.

###  VII-C Dynamic Point Cloud Analysis

Figure 9: Example of Transformed Signals in a Dynamic Point Cloud.

Spectral analysis of signals is one of the basic tools in data analysis. 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 network 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 [40]. 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 an
Euclidean threshold $\tau_{s}$ [6]. Let $\mathbf{x}_{i}$ be the 3D coordinates
of the $i$th point. We assign an edge weight between two points
$\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ as a nonzero
$A_{ij}=\exp(-\|\mathbf{x}_{i}-\mathbf{x}_{j}\|^{2}_{2}/{\sigma^{2}})$ only if
$\|\mathbf{x}_{i}-\mathbf{x}_{j}\|^{2}_{2}\leq\tau_{s}$. Next, we estimate the
spatial and temporal basis vectors of each shorter-term sequences by HOSVD in
Eq. (25). 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 short-time Fourier
transform (STFT) [41]. In Fig. 9, 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. 9, one can easily identify
different motions based on the MLN singular analysis.

Our future works shall target more interpretable analysis of these works. For
example, the physical meaning of nodes can be identified via filtered signals
under MLN highpass filter, shown in Appendix F.

##  VIII Conclusion

In this work, we present a novel tensor-based framework of multilayer network
signal processing (M-GSP) that naturally generalizes the traditional GSP to
multilayer networks. We first present the basic foundation and definitions of
M-GSP including MLN signals, signal shifting, spectral 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.

## References

  * [1] A. Sandryhaila, and J. M. F. Moura, “Discrete signal processing on graphs” IEEE Transactions on Signal Processing, vol. 61, no. 7, pp. 1644-1656, Apr. 2013.
  * [2] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega and P. Vandergheynst, “The emerging field of signal processing on graphs: extending high-dimensional data analysis to networks and other irregular domains,” in IEEE Signal Processing Magazine, vol. 30, no. 3, pp. 83-98, May 2013.
  * [3] A. Ortega, P. Frossard, J. Kovačević, J. M. F. Moura and P. Vandergheynst, “Graph signal processing: overview, challenges, and applications,” in Proceedings of the IEEE, vol. 106, no. 5, pp. 808-828, May 2018.
  * [4] A. Sandryhaila and J. M. F. Moura, ”Discrete signal processing on graphs: frequency analysis,” in IEEE Transactions on Signal Processing, vol. 62, no. 12, pp. 3042-3054, Jun., 2014.
  * [5] S. Chen, A. Sandryhaila, J. M. F. Moura and J. Kovacevic, “Signal denoising on graphs via graph filtering,” 2014 IEEE GlobalSIP, Atlanta, GA, USA, Dec. 2014, pp. 872-876.
  * [6] S. Chen, D. Tian, C. Feng, A. Vetro and J. Kovačević, “Fast resampling of three-dimensional point clouds via graphs,” in IEEE Transactions on Signal Processing, vol. 66, no. 3, pp. 666-681, Feb., 2018.
  * [7] S. Chen, A. Sandryhaila, J. M. F. Moura and J. Kovačević, “Adaptive graph filtering: multiresolution classification on graphs,” 2013 IEEE Global Conf. on Signal and Info. Processing, Austin, TX, USA, Dec. 2013, pp. 427-430.
  * [8] S. Zhang, Z. Ding and S. Cui, “Introducing hypergraph signal processing: theoretical foundation and practical applications,” in IEEE Internet of Things Journal, vol. 7, no. 1, pp. 639-660, Jan. 2020.
  * [9] S. Barbarossa and S. Sardellitti, “Topological signal processing over simplicial complexes,” in IEEE Transactions on Signal Processing, vol. 68, pp. 2992-3007, Mar. 2020.
  * [10] M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter, “Multilayer networks,” in Journal of complex networks, vol. 2, no. 3, pp. 203-271, Jul. 2014.
  * [11] M. De Domenico, A. Sole-Ribalta, E. Cozzo, M. Kivela, Y. Moreno, M. A. Porter, S. Gomez, and A. Arenas, “Mathematical formulation of multilayer networks,” Physical Review X, vol. 3, no. 4, p. 041022, 2013.
  * [12] S. Zhang, H. Zhang, H. Li and S. Cui, “Tensor-based spectral analysis of cascading failures over multilayer complex systems,” in Proc. 56th Allerton Conf. on Communication, Control, and Computing, Monticello, USA, Oct. 2018, pp. 997-1004.
  * [13] Z. Huang, C. Wang, M. Stojmenovic, and A. Nayak, “Characterization of cascading failures in interdependent cyberphysical systems,” IEEE Transactions on Computers, vol. 64, no. 8, pp. 2158-2168, Aug. 2015.
  * [14] S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, and S. Havlin, “Catastrophic cascade of failures in interdependent networks,” Nature, vol. 464, no. 7291, pp. 1025-1028, Apr. 2010.
  * [15] S. Gomez, A. Diaz-Guilera, J. Gomez-Gardenes, C. J. Perez-Vicente, Y. Moreno, and A. Arenas, “Diffusion dynamics on multiplex networks,”, Physical Review Letters, vol. 110, no. 2, p. 028701, Jan. 2013.
  * [16] P. Das and A. Ortega, “Graph-based skeleton data compression,” 2020 IEEE 22nd International Workshop on Multimedia Signal Processing (MMSP), Tampere, Finland, Sep. 2020, pp. 1-6.
  * [17] F. Grassi, A. Loukas, N. Perraudin and B. Ricaud, “A time-vertex signal processing framework: scalable processing and meaningful representations for time-series on graphs,” in IEEE Trans. Signal Processing, vol. 66, no. 3, pp. 817-829, Feb. 2018
  * [18] J. S. Stanley, E. C. Chi and G. Mishne, “Multiway graph signal processing on tensors: integrative analysis of irregular geometries,” in IEEE Signal Processing Magazine, vol. 37, no. 6, pp. 160-173, 2020.
  * [19] M. A. Porter, “What is a multilayer network,” Notices of the AMS, vol. 65, no. 11, pp.1419- 1423, Dec. 2018.
  * [20] M. De Domenico, V. Nicosia, A. Arenas, and V. Latora, “Structural reducibility of multilayer networks,” Nature Communications, vol. 6, no. 6864, pp. 1-9, Apr. 2015.
  * [21] S. Chen, R. Varma, A. Sandryhaila and J. Kovačević, “Discrete Signal Processing on Graphs: Sampling Theory,” in IEEE Transactions on Signal Processing, vol. 63, no. 24, pp. 6510-6523, Dec.15, 2015.
  * [22] A. Sandryhaila, and J. M. F. Moura, “Discrete signal processing on graphs: graph filters,” in Proceedings of 2013 IEEE ICASSP, Vancouver, Canada, May 2013, pp. 6163-6166.
  * [23] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM Review, vol. 51, no. 3, pp. 455-500, Aug. 2009.
  * [24] H. A. L. Kiers, “Towards a standardized notation and terminology in multiway analysis,” Journal of Chemometrics: A Journal of the Chemometrics Society, vol. 14, no. 3, pp. 105-122, Jan. 2000.
  * [25] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinear singular value decomposition,” SIAM Journal on Matrix Analysis and Applications, vol. 21, no. 4, pp. 1253-1278, Jan. 2000.
  * [26] A. Afshar, J. C. Ho, B. Dilkina, I. Perros, E. B. Khalil, L. Xiong, and V. Sunderam, “Cp-ortho: an orthogonal tensor factorization framework for spatio-temporal data,” Proc. 25th ACM SIGSPATIAL Intl. Conf. Advances in Geographic Info. Syst., Redondo Beach, USA, Jan. 2017, p. 67.
  * [27] I. V. Oseledets, “Tensor-train decomposition,” SIAM Journal on Scientific Computing, vol. 33, no. 5, pp. 2295-2317, Jun. 2011.
  * [28] A. Bogojeska, S. Filiposka, I. Mishkovski and L. Kocarev, ”On opinion formation and synchronization in multiplex networks,” 2013 21st TELFOR, Belgrade, Serbia, Nov. 2013, pp. 172-175.
  * [29] A. C. Kinsley, G. Rossi, M. J. Silk, and K. VanderWaal, “Multilayer and multiplex networks: an introduction to their use in veterinary epidemiology,” Frontiers in Veterinary Science, vol. 7, p. 596, Sep. 2020.
  * [30] M. De Domenico, C. Granell, M. A. Porter, and A. Arenas, “The physics of spreading processes in multilayer networks,” Nature Physics, vol. 12, no. 10, pp. 901-906, Aug. 2016.
  * [31] A. Loukas, “Graph reduction with spectral and cut guarantees,” Journal of Machine Learning Research, vol. 20, no. 116, pp. 1-42, Jun. 2019.
  * [32] Q. Shi, Y. -M. Cheung, Q. Zhao and H. Lu, “Feature extraction for incomplete data via low-rank tensor decomposition with feature regularization,” in IEEE Transactions on Neural Networks and Learning Systems, vol. 30, no. 6, pp. 1803-1817, Jun. 2019.
  * [33] M. Jouni, M. D. Mura and P. Comon, “Hyperspectral image classification using tensor cp decomposition,” IEEE International Geoscience and Remote Sensing Symposium, Yokohama, Japan, Jul. 2019, pp. 1164-1167.
  * [34] A. Sandryhaila, and J. M. F. Moura, “Discrete signal processing on graphs: graph filters,” in Proceedings of 2013 IEEE ICASSP, Vancouver, Canada, May 2013, pp. 6163-6166.
  * [35] S. Chen, F. Cerda, P. Rizzo, J. Bielak, J. H. Garrett and J. Kovačević, “Semi-supervised multiresolution classification using adaptive graph filtering with application to indirect bridge structural health monitoring,” in IEEE Trans. Signal Processing, 62(11):2879-2893, Jun. 2014.
  * [36] S. Mallat, A Wavelet Tour of Signal Processing, 3rd ed. New York, NY, USA: Academic, 2008.
  * [37] U. V. Luxburg, “A tutorial on spectral clustering.” Statistics and computing, vol .17, no. 4, pp. 395-416, 2007.
  * [38] H. Abdi, and L. J. Williams, “Principal component analysis,” Wiley Interdisciplinary Reviews: Computational Statistics, vol. 2, no. 4, pp. 433-459, 2010.
  * [39] D. Martin, C. Fowlkes, D. Tal and J. Malik, “A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics,” Proceedings Eighth IEEE ICCV, Vancouver, BC, Canada, Jul. 2001, pp. 416-423.
  * [40] CMU, “Carnegie Mellon University Graphics Lab: Motion Capture Database,” accessed on Apr. 21, 2016. [Online].
  * [41] L. Durak, and O. Arikan, “Short-time fourier transform: two fundamental properties and an optimal implementation,” IEEE Transactions on Signal Processing, vol. 51, no. 5, pp. 1231-1242, May 2003.

(a) CP Decomposition.

(b) Tucker Decomposition.

Figure 10: Diagram of Tensor Decomposition.

## Appendix

Since the initial submission is limited to 13 pages and the manuscript is
allowed to grow up in the revision stages of review, we provide additional
examples and illustrations in this Appendix. Note that the previous manuscript
should be already self-contained, and these supplementary materials are for a
better interpretation of M-GSP.

###  VIII-A Proof of Properties

####  VIII-A1 Proof of Property 1

###### Proof.

Suppose $(\lambda,\mathbf{x})$ is an eigenpair of $\mathbf{A}_{FL}$, i.e.,

$\mathbf{A}_{FL}\cdot\mathbf{x}=\lambda\mathbf{x}.$ (44)

Let $x_{N(\alpha-1)+i}=y_{M(i-1)+\alpha}$. Since

$\displaystyle A_{\alpha i\beta j}$
$\displaystyle=[A_{FL}]_{N(\alpha-1)+i,N(\beta-1)+j}$
$\displaystyle=[A_{FN}]_{M(i-1)+\alpha,M(j-1)+\beta},$ (45)

we have

$\mathbf{A}_{FN}\cdot\mathbf{y}=\lambda\mathbf{y}.$ (46)

Thus, $\lambda$ is also an eigenvalue of $\mathbf{A}_{FN}$. ∎

####  VIII-A2 Proof of Property 3

###### Proof.

Suppose that $\tilde{\mathbf{V}}_{\alpha i}$ is one factor tensor of
$\mathbf{F}$ obtained from Eq. (16).

Let $\delta[k]$ denote the Kronecker delta. Since $\mathbf{f}_{\alpha}$ forms
an orthonormal basis, then the inner product would satisfy

$<\mathbf{f}_{\alpha},\mathbf{f}_{\beta}>=\sum_{k}[f_{\alpha}]_{k}\cdot[f_{\beta}]_{k}=\delta[\alpha-\beta].$

Similarly,

$<\mathbf{e}_{i},\mathbf{e}_{j}>=\delta[i-j].$
$\displaystyle[\mathbf{F}\diamond\tilde{\mathbf{V}}_{\alpha i}]_{\beta j}$
$\displaystyle=\sum_{\sigma=1}^{M}\sum_{k=1}^{N}F_{\beta j\sigma
k}\;[\tilde{V}_{\alpha i}]_{\sigma k}$ (47)
$\displaystyle\approx\sum_{\sigma=1}^{M}\sum_{k=1}^{N}\sum_{\gamma=1}^{M}\sum_{t=1}^{N}\lambda_{\gamma
t}[f_{\gamma}]_{\beta}\;[e_{t}]_{j}\;[f_{\gamma}]_{\sigma}\;[e_{t}]_{k}\;[\tilde{V}_{\alpha
i}]_{\sigma k}.$ (48)

Then, we have

$\displaystyle\sum_{\sigma=1}^{M}\sum_{k=1}^{N}[f_{\gamma}]_{\sigma}\;[e_{t}]_{k}\;[\tilde{V}_{\alpha
i}]_{\sigma k}$
$\displaystyle=\sum_{\sigma=1}^{M}\sum_{k=1}^{N}[f_{\gamma}]_{\sigma}\;[e_{t}]_{k}\;[f_{\alpha}]_{\sigma}\;[e_{i}]_{k}$
$\displaystyle=\sum_{\sigma=1}^{M}[f_{\gamma}]_{\sigma}\;[f_{\alpha}]_{\sigma}\;\sum_{k=1}^{N}[e_{t}]_{k}\;[e_{i}]_{k}$
$\displaystyle=<\mathbf{f}_{\gamma},\mathbf{f}_{\alpha}>\cdot<\mathbf{e}_{t},\mathbf{e}_{i}>$
$\displaystyle=\delta[\gamma-\alpha]\delta[t-i]$ (49)

Thus,

$[\mathbf{F}\diamond\tilde{\mathbf{V}}_{\alpha i}]_{\beta
j}\approx\lambda_{\alpha i}[f_{\alpha}]_{\beta}[e_{i}]_{j},$ (50)

which indicates

$\mathbf{F}\diamond\tilde{\mathbf{V}}_{\alpha i}\approx\lambda_{\alpha
i}\tilde{\mathbf{V}}_{\alpha i}.$ (51)

Then, $\tilde{\mathbf{V}}_{\alpha i}$ is the approximated eigen-tensor. ∎

####  VIII-A3 Proof of Property 4

###### Proof.

Let $\tilde{\mathbf{V}}_{k}=\tilde{\mathbf{V}}_{\alpha
i}=\mathbf{f}_{\alpha}\circ\mathbf{e}_{i}$ and
$\tilde{\mathbf{V}}_{t}=\tilde{\mathbf{V}}_{\beta
j}=\mathbf{f}_{\beta}\circ\mathbf{e}_{j}$. Then, we have

$\displaystyle W_{\alpha i\beta j}$
$\displaystyle=\sum_{\theta,p}[\tilde{V}_{k}]_{\theta
p}[\tilde{V}_{t}]_{\theta p}$
$\displaystyle=\sum_{p}[e_{i}]_{p}[e_{j}]_{p}\sum_{\theta}[f_{\alpha}]_{\theta}[f_{\beta}]_{\theta}$
$\displaystyle=<\mathbf{f}_{\alpha},\mathbf{f}_{\beta}>\cdot<\mathbf{e}_{i},\mathbf{e}_{j}>,$
$\displaystyle=\delta[\alpha-\beta]\delta[i-j].$ (52)

∎

###  VIII-B Tensor Decomposition and M-GSP

####  VIII-B1 Examples and Interpretation of Tensor Decomposition

CANDECOMP/PARAFAC (CP) decomposition decomposes a tensor as a sum of the
tensor outer product of rank-one tensors [R1]. For example, a third order
tensor $\mathbf{T}\in\mathbb{R}^{I\times J\times K}$ is decomposed by CP
decomposition into

$\mathbf{T}\approx\sum_{r=1}^{R}\mathbf{a}_{r}\circ\mathbf{b}_{r}\circ\mathbf{c}_{r},$
(53)

where $\mathbf{a}_{r}\in\mathbb{R}^{I}$, $\mathbf{b}_{r}\in\mathbb{R}^{J}$,
$\mathbf{c}_{r}\in\mathbb{R}^{K}$ and a positive integer $R$ is its rank,
i.e., the smallest number of rank-one tensors in the decomposition. The
process of CP decomposition for a third-order tensor is illustrated in Fig.
10(a), and could be interpreted as the 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 [R1].
Defining a core tensor $\mathbf{S}=[S_{pqr}]\in\mathbb{R}^{P\times Q\times
R}$. The Tucker decomposition of a third-order tensor
$\mathbf{T}\in\mathbb{R}^{I\times J\times K}$ is

$\displaystyle\mathbf{T}$
$\displaystyle\approx\mathbf{S}\times_{1}\mathbf{A}\times_{2}\mathbf{B}\times_{3}\mathbf{C}$
(54)
$\displaystyle=\sum_{p=1}^{P}\sum_{q=1}^{Q}\sum_{r=1}^{R}S_{pqr}\mathbf{a}_{p}\circ\mathbf{b}_{q}\circ\mathbf{c}_{r},$
(55)

where
$\mathbf{A}=[\mathbf{a}_{1}\;\cdots\;\mathbf{a}_{p}]\in\mathbb{R}^{I\times
P}$,
$\mathbf{B}=[\mathbf{b}_{1}\;\cdots\;\mathbf{b}_{Q}]\in\mathbb{R}^{J\times
Q}$,
$\mathbf{C}=[\mathbf{c}_{1}\;\cdots\;\mathbf{c}_{R}]\in\mathbb{R}^{K\times
R}$. The diagram of the Tucker decomposition for a third-order tensor is shown
in Fig. 10(b). Note that Tucker decomposition reduces to CP decomposition if
the core tensor is limited to be super-diagonal.

####  VIII-B2 Calculation of Tensor Decomposition

To calculate the CP decomposition, approximated algorithms are implemented by
minimizing the error between the $N$th order original tensor $\mathbf{X}$ and
the factorized vectors, i.e.,

$\min_{\lambda,\mathbf{f}}||\mathbf{X}-\sum_{r=1}^{R}\lambda_{r}\mathbf{f}_{r}^{(1)}\circ\mathbf{f}_{r}^{(2)}\circ\cdots\mathbf{f}_{r}^{(N)}||^{2}_{F}.$
(56)

In addition, constraints such as orthogonality of $\mathbf{f}_{r}$ and
$||\mathbf{f}_{r}^{(i)}||^{2}=1$ can be considered during the optimization
process [R2]. Typical algorithms include alternating least squares and
projected gradient descent. Compared to CP decomposition, HOSVD is more stable
and robust. Interested readers could refer to [R3] for the implementation of
HOSVD.

TABLE III: Error of Decomposing the Representing Tensor Graph Structure | ER(0.3,0.3,6,5) | ER(0.5,0.7,11,15) | ER(0.8,0.7,6,15)
---|---|---|---
Eigen-tensor | 8.3893e-15 | 1.6001e-13 | 3.8347e-13
HOSVD | 1.011e-14 | 1.9563e-13 | 1.9056e-13
OPT-CP | 9.22e-01 | 8.82e-01 | 9.24e-01
Tucker | 9.37e-05 | 9.10e-05 | 9.40e-05

####  VIII-B3 Comparison of Different Decomposition Methods

To compare the accuracy in recovering the representing tensor using different
tensor decomposition methods we examine eigen-tensor decomposition, HOSVD,
optimal CP decomposition and Tucker decomposition in MLNs randomly generated
from the $Erd\ddot{o}s-R\dot{e}nyi$ (ER) random graph $ER(p,q,M,N)$, where $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 the decomposed and the original tensors. From the results of
Table III, we can see that the eigen-tensor decomposition and HOSVD exhibit
better overall accuracy. Generally, eigen-tensor decomposition is more
suitable for applications emphasizing joint features of layers and entities.
On the other hand, HOSVD is more efficient in separating individual features
of layers and entities. Note that, here we only used recovery accuracy to
measure different decompositions. However, different decompositions may have
different performance advantages when capturing different data features that
can be measured with different metrics.

###  VIII-C Spectral Space of Directed MLN

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 spectral space of a directed multilayer network,
we also present two ways to compute: 1) Flattening analysis; and 2) Tensor
factorization:

  * •

Flattening analysis: Similar to the representing tensor of undirected graphs,
we flatten the representing tensor as a second-order supra-matrix, and define
spectral space as the reshape of the eigenvectors of the supra-matrix. The
flattened matrix $\mathbf{A}_{FX}$ (or $\mathbf{A}_{FN}$, $\mathbf{A}_{FL}$)
can be decomposed as

$\mathbf{A}_{FX}\approx\mathbf{E}\Sigma\mathbf{E}^{-1},$ (57)

where $\mathbf{E}\in\mathbb{R}^{MN\times MN}$ is the matrix of eigenvectors
and $\Sigma=diag(\lambda_{i})$ is a diagonal matrix of eigenvalues. Then, we
can reshape the eigenvectors, i.e., each column of $\mathbf{E}$ as
$\mathbf{V}_{k}\in\mathbb{R}^{M\times N}$, and reshape each row of
$\mathbf{E}^{-1}$ as $\mathbf{U}_{k}\in\mathbb{R}^{M\times N}$. Consequently,
the original tensor can decomposed into

$\mathbf{A}\approx\sum_{k=1}^{MN}\lambda_{k}\mathbf{V}_{k}\circ\mathbf{U}_{k}.$
(58)

  * •

Tensor Factorization: We can also compute the spectrum from the tensor
factorization based on CP-decomposition

$\displaystyle\mathbf{A}$
$\displaystyle\approx\sum_{k=1}^{R}\lambda_{k}\mathbf{a}_{k}\circ\mathbf{b}_{k}\circ\mathbf{c}_{k}\circ\mathbf{d}_{k}$
(59)
$\displaystyle=\sum_{k=1}^{R}\lambda_{k}\mathbf{V}_{k}\circ\mathbf{U}_{k}.$
(60)

where $R$ is the rank of tensor,
$\mathbf{a}_{k},\mathbf{c}_{k}\in\mathbb{R}^{M}$ characterize layers,
$\mathbf{b}_{k},\mathbf{d}_{k}\in\mathbb{R}^{N}$ characterize entities, and
$\mathbf{V}_{k}=\mathbf{a}_{k}\circ\mathbf{b}_{k},\mathbf{U}_{k}=\mathbf{c}_{k}\circ\mathbf{d}_{k}\in\mathbb{R}^{M\times
N}$ characterize the joint features. Since there are $MN$ nodes, it is clear
that $R\leq MN$. Note that, for a single layer, Eq. (59) reduces to

$\displaystyle\mathbf{A}\approx\sum_{k=1}^{N}\lambda_{k}\mathbf{v}_{k}\circ\mathbf{u}_{k}.$
(61)

Moreover, if $\mathbf{V}=(\mathbf{v}_{k})$ and
$\mathbf{U}=(\mathbf{u}_{k}^{\mathrm{T}})=\mathbf{V}^{-1}$ are orthogonal
bases, Eq. (61) is in a consistent form of the eigendecomposition in a single-
layer normal graph. In addition, Eq. (16) is also a special case of Eq. (59)
if the multilayer network 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.

###  VIII-D Interpretation of Total Variation in DSP

In DSP, the discrete Fourier transform (DFT) of a sequence $s_{n}$ is given by
$\hat{s}_{k}=\sum_{n=0}^{N-1}s_{n}e^{-j\frac{2\pi kn}{N}}$ and the frequency
is defined as $\nu_{n}=\frac{n}{N}$, $n=0,1,\cdots,N-1$. We can easily
summarize the following conclusions:

  * •

$\nu_{n}:\;1<n<\frac{N}{2}-1$ corresponds to a continuous time signal
frequency ${n\over N}f_{s}$;

  * •

$\nu_{n}:\;\frac{N}{2}+1<n<N-1$ corresponds to a continuous time signal
frequency $-(1-{n\over N})f_{s}$;

  * •

$\nu_{\frac{N}{2}}$ corresponds to $f_{s}/2$;

  * •

$n=0$ corresponds to frequency 0.

Here, $f_{s}$ is the critical sampling frequency. In traditional DFT, we
generate the Fourier transform
$\hat{f}(\omega)=\int_{-\infty}^{\infty}f(x)e^{-2\pi jx\omega}dx$ at each
discrete frequency $\frac{n}{N}f_{s}$,
$n=-\frac{N}{2}+1,-\frac{N}{2}+2,\cdots,\frac{N}{2}-1,\frac{N}{2}$. The
highest and lowest frequencies correspond to $n=N/2$ and $n=0$, respectively.
Note that $n$ varies from $-\frac{N}{2}+1$ to $\frac{N}{2}$ here. Since
$e^{-j2\pi k\frac{n}{N}}=e^{-j2\pi k\frac{n+N}{N}}$, we can let $n$ vary from
$0$ to $N-1$ and cover the complete period. Now, $n$ varies in exact
correspondence to $\nu_{n}$, and the aforementioned conclusions are drawn. The
highest frequency occurs at $n=\frac{N}{2}$.

The total variation in DSP is defined as the differences among the signals
over time [R4], i.e.,

$\displaystyle\mathbf{TV(s)}$
$\displaystyle=\sum_{n=0}^{N-1}|s_{n}-s_{(n-1\mod N)}|$ (62a)
$\displaystyle=\mathbf{||s-C_{N}s||}_{1},$ (62b)

where

$\displaystyle{\mathbf{C_{N}}=\begin{bmatrix}0&0&\cdots&0&1\\\
1&0&\cdots&0&0\\\ \vdots&\ddots&\ddots&\ddots&\vdots\\\ 0&0&\ddots&0&0\\\
0&0&\cdots&1&0\end{bmatrix}.}$ (63)

When we perform the eigen-decomposition of $\mathbf{C_{N}}$, we see that the
eigenvalues are $\lambda_{n}=e^{-j\frac{2\pi n}{N}}$ with eigenvector
$\mathbf{f}_{n}$, $0\leq n\leq N-1$. More specifically, the total variation of
the frequency component $\mathbf{f}_{n}$ is calculated as

$\mathbf{TV(f}_{n})=|1-e^{-j\frac{2\pi n}{N}}|,$ (64)

which increases with $n$ for $n\leq\frac{N}{2}$ before decreasing with $n$ for
$\frac{N}{2}<n\leq N-1$.

(a) Original Image.

(b) K-Means.

(c) GSP.

(d) MLN-SVD.

(e) MLN-FZ.

Figure 11: Additional Visualization Results. TABLE IV: Results of Image Segmentation in Image BSD300 and BSD500 | BSD300(N=100) | BSD300(N=300) | BSD300(N=100, Coarse) | BSD300(N=300, Coarse) | BSD300(N=900,Coarse) | BSD500(Coarse)
---|---|---|---|---|---|---
GSP | 0.1237 | 0.1149 | 0.3225 | 0.3087 | 0.3067 | 0.3554
K-MEANS | 0.1293 | 0.1252 | 0.3044 | 0.3105 | 0.3124 | 0.3154
MLN-SVD |  0.1326 |  0.1366 |  0.3344 |  0.3394 |  0.3335 |  0.3743
MLN-CP | 0.1321 | 0.1293 | 0.3195 | 0.3256 | 0.3243 | 0.3641
IIC |  | 0.2071
GS | 0.3658
BP | 0.3239
DFC | 0.3739

###  VIII-E Filter Design Optimization

Filters designed through optimization based on geometric information play an
important role in semi-supervised signal processing works, such as
classification, denoising, and prediction [R5-R7]. We can also design filters
based on similar optimization approaches. Let an MLN optimization-based filter
be denoted by $\mathbf{s}^{\prime}=f(\mathbf{s},\mathbf{F})$ where
$\mathbf{s}$ is the input signal. The filter design for $f(\cdot)$ can be
posed as an optimization formulation:

$\min_{f}\quad F_{1}(\mathbf{s},\mathbf{s}^{\prime})+\sigma
F_{2}(\mathbf{s}^{\prime},\mathbf{F}).$ (65)

Here, $F_{1}$ denotes a chosen empirical loss, such as the least square loss
or the hinge loss, whereas $F_{2}$ is a regularization loss, such as Laplacian
regularization and total variation. In addition to optimization when designing
filter $\mathbf{f}(\cdot)$, some applications focus on optimization of the
filtered signal $\mathbf{s}^{\prime}$ directly, which manifests a similar form
as the traditional graph-based semi-supervised learning approaches [R8].

###  VIII-F Additional Results

####  VIII-F1 Results of Image Segmentation

Here, we provide additional visualization results of the BCCD datasets. In
addition, numerical result in BSD500 are also provided, compared to learning-
based methods including invariant information clustering (IIC) [R9], graph-
based segmentation (GS) [R10], back propagation (BP) [R11] and differentiable
feature clustering (DFC) [R12] in Table IV Shown as Fig. 11, MLN-based methods
show a more robust result than k-means and GSP-based methods. For the coarse
scenarios with fewer clusters, geometry-based methods generally work better.
MLN 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 [R12].

####  VIII-F2 Dynamic Point Cloud Analysis via MLN Highpass Filter

To understand the motions in the dynamic point clouds, we can also apply the
entity-wise MLN highpass filters described in Section V-B to capture the key
details of human bodies. More specifically, we pick the first 140 frames in
‘walking’ and define the norm of three coordinates as signals. We select 5
joints (entities) of the body in each temporal frame (layers) shown as Fig.
12. From the results, entity 1 and entity 2 display periodic patterns which
are related to the leg movement. Entity 3 almost show no changes since the
head does not move compared to the main body. Entity 4 and entity 5 display
more irregular patterns since hands do not directly identify the pose of
‘walking’. To summarize, the MLN highpass filter can efficiently capture the
key information in the motions and identify the meaning of nodes (entities).

Figure 12: Example of Filtered Signals in a Dynamic Point Cloud: entity 1 and
entity 2 are legs; entity 3 is head; entity 4 and entity 5 are hands.

## References

  * [R1] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM Review, vol. 51, no. 3, pp. 455-500, Aug. 2009.
  * [R2] A. Afshar, J. C. Ho, B. Dilkina, I. Perros, E. B. Khalil, L. Xiong, and V. Sunderam, “Cp-ortho: an orthogonal tensor factorization framework for spatio-temporal data,” in Proceedings of the 25th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems, Redondo Beach, CA, USA, Jan. 2017, p. 67.
  * [R3] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinear singular value decomposition,” SIAM Journal on Matrix Analysis and Applications, vol. 21, no. 4, pp. 1253-1278, Jan. 2000.
  * [R4] S. Mallat, A Wavelet Tour of Signal Processing, 3rd ed. New York, NY, USA: Academic, 2008.
  * [R5] S. Chen, F. Cerda, P. Rizzo, J. Bielak, J. H. Garrett and J. Kovačević, “Semi-supervised multiresolution classification using adaptive graph filtering with application to indirect bridge structural health monitoring,” in IEEE Transactions on Signal Processing, vol. 62, no. 11, pp. 2879-2893, Jun. 2014.
  * [R6] K. Avrachenkov, A. Mishenin, P. Gonçalves, and M. Sokol, “Generalized optimization framework for graph-based semi-supervised learning,” in Proceedings of the 2012 SIAM International Conference on Data Mining, Anaheim, CA, USA, Apr. 2012, pp. 966-974.
  * [R7] J. Kunegis, S. Schmidt, A. Lommatzsch, J. Lerner, E. W. De Luca, and S. Albayrak, “Spectral analysis of signed graphs for clustering, prediction and visualization,” in Proceedings of the 2010 SIAM International Conference on Data Mining, Columbus, OH, USA Apr. 2010, pp. 559-570.
  * [R8] D. Zhou, J. Huang, and B. Schölkopf, “Learning with hypergraphs: clustering, classification, and embedding,” Advances in Neural information Processing Systems, Vancouver, Canada, Dec. 2006, pp. 1601-1608.
  * [R9] X. Ji, J. F. Henriques, and A. Vedaldi, “Invariant information clustering for unsupervised image classification and segmentation,” in Proceedings of the IEEE International Conference on Computer Vision, Seoul, Korea, Nov. 2019, pp. 9865– 9874\.
  * [R10] P. F. Felzenszwalb and D. P. Huttenlocher, “Efficient graph-based image segmentation,” International Journal of Computer Vision, vol. 59, no. 2, pp. 167–181, 2004.
  * [R11] A. Kanezaki, “Unsupervised image segmentation by backpropagation,” in Proceedings of IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP), Calgary, AB, Canada, Apr. 2018, pp. 1543-1547.
  * [R12] W. Kim, A. Kanezaki and M. Tanaka, “Unsupervised learning of image segmentation based on differentiable feature clustering,” in IEEE Transactions on Image Processing, vol. 29, pp. 8055-8068, Jul. 2020.