\section{Introduction}\label{intro}
\IEEEPARstart{G}{eometric} signal processing tools have found broad applications in data analysis to uncover obscure or hidden structures from complex datasets \cite{c1}.
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
\cite{c3}. For a signal with $N$ samples, a graph of $N$ nodes can be formed
to model their underlying interactions \cite{c2}. In GSP, a graph Fourier space is also
defined from the spectral space of the representing matrix (adjacency/Laplacian)
for signal processing tasks \cite{c4}, such as denoising \cite{c5}, resampling \cite{c6}, and classification \cite{c7}. Generalization of the more traditional GSP
includes signal processing over hypergraphs \cite{c8} and simplicial complexes \cite{c9},
{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 \cite{c11}. For example, in a video
dataset modeled by spatial temporal graph shown in Fig. \ref{ex1},
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 \cite{c10}.

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. \ref{ex2}: the power grid and the computation network. These two layers have different physical connectivity and rules \cite{c12}.
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 \cite{c13}.
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 \cite{c14}.

\begin{figure*}[htbp]
	\centering
	\subfigure[]{
		\label{ex1}

		\includegraphics[height=2.5cm]{video.png}}
	\hspace{3cm}
	\subfigure[]{
		\label{ex2}
		\includegraphics[height=2.5cm]{CPS.png}}
	\caption{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.}
	\label{ex_mln}
	\vspace{-3mm}
\end{figure*}

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 \cite{d1}, 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 \cite{c16}, 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 \cite{c19}, 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 \cite{c18}. 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 \cite{d2}, 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 \cite{c20}.
However, such matrix-based representations mainly focus on the intralayer connections and lack representation for interlayer interactions \cite{c20}.
A more natural and general way may start with tensor representation
\cite{c10}, 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 {\em intralayer} and {\em 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:
\begin{itemize}
	\item Leveraging tensor representation of MLN, we introduce the definitions of signals and shifting over MLN.
	\item 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.
	\item We also present  fundamentals of filter design in M-GSP, and suggest several practical applications based on the proposed framework.
\end{itemize}

We organize the technical contents as follows.
Section \ref{prelim} first summarizes the
preliminaries of traditional GSP and tensor analysis.
We then introduce the fundamentals
of M-GSP and frequency analysis in Section \ref{funda} and Section \ref{mln_spec}, respectively. We next present MLN filter design in Section \ref{fter}. We provide the physical insights and
spectrum interpretation of M-GSP concepts in Section \ref{discus}. Within the proposed M-GSP framework, we present several example applications and demonstrate its effectiveness in Section \ref{app}, before summarizing our conclusions in Section \ref{con}.

\section{Preliminaries} \label{prelim}
\subsection{Overview of Graph Signal Processing}
Signal processing on graphs \cite{c1,c2,c3} 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}'=\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
\begin{equation}\label{GFT}
\hat{\mathbf{s}}=\mathbf{V}^{-1}\mathbf{s},
\end{equation}
whereas the inverse GFT is given by
$\mathbf{s}=\mathbf{V}\hat{\mathbf{s}}$.

From definitions of GFT, other concepts, such as
sampling theory \cite{c23}, filter design \cite{c24}, and frequency analysis \cite{c4}
can be developed for signal processing and data analysis tasks.

\subsection{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 \cite{c25}.
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 \cite{c25}.
\subsubsection{Super-diagonal Tensor}
An $N$th-order tensor $\mathbf{A}\in\mathbb{R}^{I_1\times I_2\cdots\times I_N}$ is \textit{super-diagonal} if its entries $A_{i_1i_2\cdots i_N}\neq 0$ only for $i_1=i_2=\cdots=i_N$.

\subsubsection{Symmetric Tensor}
A tensor is \textit{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 \textit{super-symmetric} if
$A_{ijk}=A_{jik}=A_{kij}=A_{kji}=A_{jik}=A_{jki}$.
In addition, tensors can be \textit{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 \textit{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$.

\subsubsection{Tensor Outer Product}
The \textit{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
\begin{equation}
\mathbf{W}=\mathbf{U} \circ \mathbf{V}.
\end{equation}
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.

\subsubsection{n-mode Product}
The \textit{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
\begin{equation}
\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}.
\end{equation}
 Each element in $\mathbf{W}$ is given by
$W_{i_1 i_2 \cdots i_{n-1} j i_{n+1} \cdots i_P}=\sum_{i_n=1}^{I_n}U_{i_1\cdots i_P}V_{j i_n}.$
Note that the $n$-mode product is a different operation from matrix product.

\subsubsection{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
\begin{equation}\label{contract}
	\mathbf{y}=\mathbf{A}\diamond \mathbf{x}\in \mathbb{R}^{M\times N},
\end{equation}
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
\begin{equation}\label{oodot}
	\mathbf{W}=\mathbf{U}\odot \mathbf{V}\in\mathbb{R}^{M\times
		N\times M \times N},
\end{equation}
whose entries are $W_{\alpha i \epsilon p}=\sum_{\beta j} U_{\alpha i \beta j} V_{\beta j \epsilon p}$.
\subsubsection{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 \cite{c25,c26}.
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 \cite{c25}.
Other typical decompositions include Higher-Order SVD (HOSVD) \cite{c27}, orthogonal CP-decomposition
\cite{c28}, and Tensor-Train decomposition \cite{c29}.
Interested readers are referred to the tutorial \cite{c25} for more
details. Additional examples and illustrations of tensor decomposition in M-GSP are also provided in the \textbf{Appendix B} because of page limit.

\section{Fundamentals of M-GSP} \label{funda}
In this section, we introduce the basic definitions in the proposed M-GSP framework.
\begin{figure}[t]
	\centering
	\subfigure[]{
		\label{mln1}

		\includegraphics[height=3cm]{MLN.png}}
	\hspace{2cm}
	\subfigure[]{
		\label{mlp1}
		\includegraphics[height=3cm]{MLP.png}}
	\caption{Example of multilayer networks: (a) A three-layer interconnected network; (b) a three-layer multiplex network.}
	\label{ex_mln1}
	\vspace{-5mm}
\end{figure}
\subsection{Multilayer Network}
Before introducing the foundations of M-GSP, we first provide  definitions of multilayer networks (MLN) \cite{c25}.
\begin{definition}[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.
\end{definition}
Note that, we mainly focus on the layer-disjoint multilayer network \cite{c11} 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 \cite{d3}.

In multilayer networks, edges connect nodes in the same layer (intralayer edges) or nodes of different layers (interlayer edges) \cite{c30}. There are two main types of multilayer networks: \textit{multiplex network} and \textit{interconnected network} \cite{c31}. In a \textit{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 \cite{c31}. In the \textit{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. \ref{ex_mln1}, where different colors represent different layers,  solid lines represent intralayer connections, and  dash lines indicate interlayer connections.

\subsection{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 \cite{c11}.

\subsubsection{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 \cite{c10}.
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
\begin{equation}\label{tensor_construction}
\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},
\end{equation}
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.
\begin{definition}[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
	\begin{equation}
	\mathbf{A}=(A_{\alpha i \beta j}), \quad 1\leq \alpha,\beta\leq M, 1\leq i,j\leq N.
	\end{equation}
	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$.
\end{definition}

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 \ref{mln_spec}.

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.

\begin{definition}[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.
\end{definition}

The Laplacian tensor can be useful to analyze propagation processes
such as
diffusion or random walk \cite{c10}. Both adjacency and Laplacian
tensors are important algebraic representations of the MLN depending on datasets and user objectives.

\subsubsection{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 \cite{c19} and MWGSP \cite{c18}, 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 \cite{c12} and set the augmented signals as zero; and 2) Aggregate several nodes into super-nodes for
layers with $|l_i|>N$ \cite{c33} 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. (\ref{tensor_construction}), 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 \ref{exn_la}.

\subsection{Flattening and Analysis} \label{mln_flat}
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:
\begin{itemize}
	\item Layer-wise Flattening: The representing tensor $\mathbf{F}$ can be flattened into $\mathbf{F}_{FL} \in \mathbb{R}^{MN\times MN}$ with
	each element
	\begin{equation}\label{lwise}
	{[F_{FL}]}_{{N(\alpha-1)+i, N(\beta-1)+j}}=F_{\alpha i \beta j}.
	\end{equation}
	\item Entity-wise Flattening:  The representing tensor $\mathbf{F}$ can be flattened into $\mathbf{F}_{FN} \in \mathbb{R}^{NM\times NM}$ with each element
	\begin{equation}\label{ewise}
	{[F_{FN}]}_{{M(i-1)+\alpha, M(j-1)+\beta}}=F_{\alpha i \beta j}.
	\end{equation}
\end{itemize}

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. \ref{flat}. 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
{\em 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
{\em 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 \cite{c12}, the approximations of spectral radius are derived based on different structures of these two flattened matrix.

\begin{figure}[t]
	\centering
	\includegraphics[width=3.5in]{flattening.png}
	\caption{Example of Multilayer Network Flattening.}
	\vspace{-2mm}
	\label{flat}
	\vspace{-4mm}
\end{figure}

\subsection{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 \cite{c2},  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.

\begin{definition} [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
	\begin{equation}
	\mathbf{s}=(s_{\alpha i})\in\mathbb{R}^{M\times
	N}, \quad 1\leq \alpha \leq M, 1\leq i\leq N,
	\end{equation}
	where the entry $s_{\alpha i}$ is the signal sample in the projected node of entity $i$ in layer $\alpha$.
\end{definition}

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}$:
\begin{itemize}
	\item Layer-wise flattening: $\mathbf{s}_L\in \mathbb{R}^{MN}$ whose entries are calculated as
	$[s_L]_{N(\alpha-1)+i}=s_{\alpha i}$.

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

\end{itemize}

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.

\begin{definition}[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.,
	\begin{equation} \label{shift_mln}
	\mathbf{s}'=\mathbf{F}\diamond \mathbf{s}\in \mathbb{R}^{M\times
	N},
	\end{equation}
	where $\diamond$ is the contraction between $\mathbf{F}$ and $\mathbf{s}$ defined in Eq. (\ref{contract}).
\end{definition}
The elements in the shifted signal $\mathbf{s}'$ are calculated as
\begin{equation} \label{diffuse}
s'_{\alpha i}=\sum_{\beta=1}^M\sum_{j=1}^N F_{\alpha i \beta j}s_{\beta j}.
\end{equation}

From Eq. (\ref{diffuse}), 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. (\ref{diffuse}) can be written as
\begin{equation}
	s'_{\alpha i}=\sum_{\beta=1}^M\sum_{j=1}^N A_{\alpha i \beta j} (s_{\alpha i}-s_{\beta j}),
\end{equation}
which is the weighted average of difference with neighbors.

\section{Multilayer Network Spectral Space} \label{mln_spec}
In traditional GSP,  graph spectral space is defined according to the eigenspace of the representing matrix \cite{c2}. 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 \textbf{Appendix C} and leave the frequency analysis in the future works. Meanwhile, all the proofs of the properties in this part are listed in \textbf{Appendix A}.

\subsection{Joint Spectral Analysis in M-GSP} \label{JSA}
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 \cite{c10} 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
\begin{align}\label{ted}
\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},
\end{align}
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.
\begin{definition}[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.
\end{definition}

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. (\ref{GFT}). 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.
\begin{definition}[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.,
	\begin{equation}\label{M-GFT}
	\hat{\mathbf{s}}=\mathbf{U}_\mathcal{F} \diamond \mathbf{s}.
	\end{equation}
\end{definition}

Now, we show how to obtain the eigen-tensors. Implementing the flattening
analysis, we have the following properties.

\begin{property}
	The two types of flattened tensor in Eq. (\ref{lwise}) and Eq. (\ref{ewise}) lead to the same eigenvalues.
\end{property}
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.

\begin{property}\label{eig}
	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}$.

\end{property}
The {Property \ref{eig}} 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.

\subsection{Order-wise Spectral Analysis in M-GSP}\label{OAA}

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 \cite{c28} as follows:
\begin{align}\label{decompose1}
\mathbf{F}&\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\\
&=\sum_{\alpha=1}^M\sum_{i=1}^N \lambda_{\alpha i} \tilde{\mathbf{V}}_{\alpha i} \circ\tilde{\mathbf{V}}_{\alpha i},
\end{align}
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 \cite{d5} and tensor-based PCA analysis \cite{d4}. A detailed discussion
of tensor decomposition and its implementation in M-GSP are provided in \textbf{Appendix B}.
In Eq. (\ref{decompose1}), $\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 \ref{dep}.

Note that, if there is only one layer in the multilayer network, Eq. (\ref{decompose1}) 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. (\ref{decompose1}), the order-wise MLN spectrum is defined as follows.
\begin{definition}[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.
\end{definition}

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.
\begin{property}
	The factor tensor $\tilde{\mathbf{V}}_{\alpha i}$ of the representing tensor $\mathbf{F}$ is the approximated eigen-tensor of $\mathbf{F}$.
\end{property}

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.

\begin{property}\label{ortho}
	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.
\end{property}
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. (\ref{M-GFT}), the joint transform can be defined as
\begin{equation}\label{O-GFT}
\hat{\mathbf{s}}=\tilde{\mathbf{U}}_\mathcal{F} \diamond \mathbf{s}.
\end{equation}
Note that each element of $\hat{\mathbf{s}}$ in Eq. (\ref{O-GFT})
can be calculated as
\begin{align}
\hat{s}_{\alpha i}&=\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}\\
&=\sum_{\beta,j} [f_\alpha]_\beta \cdot [e_i]_j\cdot s_{\beta j}.
\end{align}
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}}'=\mathbf{E}_f^{\mathrm{T}}\mathbf{s}\mathbf{E}_e$,
with each element $\tilde{s}_\alpha'=\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}}'=\mathbf{E}_f^{\mathrm{T}}\mathbf{s}\mathbf{E}_e$.
Then, we have the following definition of M-GFT based on order-wise spectrum.

\begin{definition} [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
	\begin{equation}
	\hat{\mathbf{s}}_L=\mathbf{E}_f^{\mathrm{T}}\mathbf{s}\in\mathbb{R}^{M\times N},
	\end{equation}
	and the entity-wise M-GFT can be defined as
	\begin{equation}
	\hat{\mathbf{s}}_N=\mathbf{s}\mathbf{E}_e\in\mathbb{R}^{M\times N}.
	\end{equation}
	The joint M-GFT based on order-wise spectrum is defined as
	\begin{equation} \label{ffff}
	\hat{\mathbf{s}}=\mathbf{E}_f^{\mathrm{T}}\mathbf{s}\mathbf{E}_e\in\mathbb{R}^{M\times N}.
	\end{equation}
\end{definition}
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. (\ref{GFT}).

In addition, since $\mathbf{f}_\alpha$ and $\mathbf{e}_i$ are orthonormal basis of undirected MLN, the inverse M-GFT can be calculated as
\begin{equation}\label{igft}
\mathbf{s}'=\mathbf{E}_f\hat{\mathbf{s}}\mathbf{E}_e^{\mathrm{T}}.
\end{equation}
Different from joint MLN Fourier space in Section \ref{JSA}, 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.

\subsection{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) \cite{c27} 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
\begin{equation}\label{decomposeS}
\mathbf{F}= \mathbf{S}\times_1 \mathbf{U}^{(1)}\times_2 \mathbf{U}^{(2)}\times_3 \mathbf{U}^{(3)}\times_4 \mathbf{U}^{(4)},
\end{equation}
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
\begin{itemize}
	\item $<\mathbf{S}_{i_n=\alpha},\mathbf{S}_{i_n=\beta}>=0$ where $\alpha\neq\beta$.
	\item $||\mathbf{S}_{i_n=1}||\geq||\mathbf{S}_{i_n=2}||\geq\cdots\geq ||\mathbf{S}_{i_n=I_n}||\geq 0$.
\end{itemize}

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
\begin{equation}\label{decom}
(\lambda_{\alpha i}, \hat{\mathbf{V}}_{\alpha i})=(\gamma_\alpha\cdot\sigma_i, \mathbf{f}_\alpha \cdot \mathbf{e}_i).
\end{equation}
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.
\begin{definition}[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. (\ref{decom}). The singular vectors $\{\mathbf{f}_1,\cdots,\mathbf{f}_M\}$ and $\{\mathbf{e}_1,\cdots,\mathbf{e}_N\}$ in Eq. (\ref{decomposeS}) characterize layers and entities, respectively.
\end{definition}

Similar to order-wise spectral analysis in Section \ref{OAA},
we can define the MLN singular tensor transform (M-GST) based on the singular tensors as follows.
\begin{definition}[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. (\ref{decomposeS}), 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.,
	\begin{equation} \label{ssss}
	\check{\mathbf{s}}=\mathbf{U}_s \diamond \mathbf{s}.
	\end{equation}
	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
	\begin{equation}
	\check{\mathbf{s}}_L=\mathbf{W}_f^{\mathrm{T}}\mathbf{s}\in\mathbb{R}^M\times \mathbb{R}^N,
	\end{equation}
	and the entity-wise M-GST can be defined as
	\begin{equation}
	\check{\mathbf{s}}_N=\mathbf{s}\mathbf{W}_e\in\mathbb{R}^M\times \mathbb{R}^N.
	\end{equation}
\end{definition}
\noindent
Inverse M-GST can be defined similarly as in Eq. (\ref{igft}) with unitary $\mathbf{W}_e$ and $\mathbf{W}_f$.

Compared to the eigen-tensors in Eq. (\ref{ted}), 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 \ref{app}, together with additional discussions in \textbf{Appendix B}. The intuition of applying HOSVD in MLN analysis and its correlations to GSP are also provided in Section \ref{i_hosvd}.

\subsection{Spectrum Ranking in the Multilayer Network} \label{ex_fre}
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 \cite{c2}.
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:
\begin{align}\label{tv}
TV(\mathbf{V}_{\alpha i})&=||\mathbf{V}_{\alpha i}-\frac{1}{|\lambda|_{max}} \mathbf{F}\diamond \mathbf{V}_{\alpha i}||_1\\
&=|1-\frac{\lambda}{|\lambda|_{max}}\mathbf|\cdot||{V}_{\alpha i}||_1,
\end{align}
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 \cite{c2}. 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 \ref{fre_int}.

\section{Filter Design} \label{fter}
In this section, we introduce an MLN filter design together with its properties based on signal shifting.

\subsection{Polynomial Filter Design}
Polynomial filters are  basic filters in
GSP \cite{c7,d6}. In M-GSP,
first-order filtering consists of
basic signal filtering, i.e.,
$\mathbf{s}'=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.,
\begin{align}
\mathbf{s}''=f_2(\mathbf{s})=\mathbf{F}\diamond(\mathbf{F}\diamond \mathbf{s}),
\end{align}
whose entries $s_{\alpha i}''$ are calculated as
\begin{align}
s_{\alpha i}''&=\sum_{\beta=1}^{M}\sum_{j=1}^N F_{\alpha i \beta j} s_{\beta j}'
=\sum_{\beta ,j} F_{\alpha i \beta j}\sum_{\epsilon ,p} F_{\beta j \epsilon p}s_{\epsilon p}\\
&=\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},
\end{align}
where $\odot$ is the contraction defined in Eq. (\ref{oodot}).

Let $\mathbf{F}^{[2]}=\mathbf{F}\odot \mathbf{F}$. From Eq. (\ref{ted}),
we have:
\begin{align}
F^{[2]}_{\alpha i\beta j}&=\sum_{\theta,p} F_{\alpha i \theta p} F_{\theta p \beta j}\\
&=\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})\nonumber\\
&=\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})\nonumber\\
&=\sum_k \lambda_k^2[ {V}_k]_{\alpha i}[ {V}_k]_{\beta j}.
\end{align}
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.
\begin{property}
	The $\tau$-th order basic shifting filter $f_\tau(\mathbf{s})$ can be calculated as
	\begin{align}\label{polyf}
	f_\tau(\mathbf{s})&=\mathbf{F}^{[\tau]}\diamond\mathbf{s}
	=(\sum_{k=1}^{MN}\lambda_k^{\tau}\mathbf{V}_k\circ\mathbf{V}_k)\diamond\mathbf{s}.
	\end{align}
\end{property}
This property is the M-GSP counterpart to the traditional
linear system interpretation
that complex
exponential signals are eigenfunction of linear systems \cite{c2}, and provides a quicker implementation of higher-order shifting.
With the $k$-order polynomial term, the adaptive polynomial filter is defined as
\begin{equation}
h(\mathbf{s})=\sum_k \alpha_k\mathbf{F}^{[k]}\diamond\mathbf{s},
\end{equation}
where $\{\alpha_k\}$ are parameters to be estimated from data.

Adaptive polynomial filter is useful in semi-supervised classification \cite{d7} and exploits
underlying geometric topologies. We will illustrate further and provide application examples based on MLN polynomial filtering in Section \ref{app}.

\subsection{Spectral Filter Design}\label{fefi}
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 \cite{c2}, 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}'=[\hat{s}_1,\cdots, \hat{s}_{N-k}, 0, \cdots, 0]$, a high-pass filter can be designed as
\begin{align}
	\mathbf{s}'&=\mathbf{V}_\mathcal{G}\hat{s}'
	=\mathbf{V}_\mathcal{G}\Sigma_k\mathbf{V}_\mathcal{G}^\mathrm{T}\mathbf{s}
\end{align}
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. (\ref{ffff}) and Eq. (\ref{ssss}),  spectral filtering of $\mathbf{s}$ is defined as
\begin{align}
	&\mathbf{s}'=\nonumber\\
&\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}}
\end{align}
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 \ref{discus} and Section \ref{app}.

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 \textbf{Appendix E} for a short discussion.

\section{Discussion and Interpretative Insights} \label{discus}
\subsection{Interpretation of M-GSP Frequency} \label{fre_int}

\subsubsection{Interpretation of Graph Frequency}\label{gfe}
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 \cite{c34}.
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 \cite{c2,c8} or \textbf{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. \ref{cirg}. Suppose that its adjacency matrix is written as
\begin{align}
{
	\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]}
\end{align}
Then, the shifted signal over the cyclic graph is calculated as $\mathbf{s}'=\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 \cite{c2}.

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.

\begin{figure}[t]
	\centering
	\includegraphics[width=1.5in]{CIRC.jpg}
	\vspace{-2mm}
	\caption{Example of Cyclic Graph.}
	\label{cirg}
		\vspace{-4mm}
\end{figure}

\subsubsection{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. (\ref{tv}). 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.
\begin{figure}[t]
	\centering
	\includegraphics[width=1.8in]{differece.jpg}
	\vspace{-2mm}
	\caption{Example of MLN Frequencies.}
	\label{diff}
		\vspace{-4mm}
\end{figure}
\begin{table}[t]
	\centering
	\caption{Total Difference of all Nodes}
	\vspace{-2mm}
	\begin{tabular}{|l|l|l|l|}
		\hline
		Eigenvalue       & 0.5726 & 2.7500 & 5.7259 \\ \hline
		Total Differences & 3.7439 & 2.6846 & 0.8293 \\ \hline
	\end{tabular}
	\label{t_diff}
	\vspace{-4mm}
\end{table}

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 \cite{c2}.
We next evaluate the differences between filtered signal and original signal by treating different eigen-tensor as signals. From  results of Fig. \ref{diff},
shifted signal sample of each node changes more than its original samples. Also from the results given n Table \ref{t_diff}, 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.

\subsubsection{Interpretation of MLN Singular Tensors}\label{i_hosvd}
As discussed in Section \ref{gfe}, 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 \cite{d11} and PCA analysis \cite{d12}. 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.

\vspace{-3mm}
\subsection{Interpretation of Entities and Layers} \label{exn_la}
To gain better physical insight of entities and layers,
we discuss two categories of datasets:
\begin{itemize}
	\item 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.

	\item 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. \ref{le},
	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.
\end{itemize}

\begin{figure}[t]
	\centering
	\includegraphics[width=2in]{le.png}
	\vspace{-1mm}
	\caption{Example of Different Entities.}
	\vspace{-5mm}
	\label{le}
\end{figure}

\vspace{-3mm}
\subsection{Distinctions from Existing GSP Works} \label{dep}
\subsubsection{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 \ref{JSA}, 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.

\subsubsection{Joint Time-Vertex Fourier Analysis}
In \cite{c19}, a joint time-vertex Fourier transform (JFT) is defined by implementing GFT and DFT consecutively.
As discussed in Section \ref{gfe}, 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.

\subsubsection{Multi-way Graph Signal Processing}
In \cite{c18}, 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.
\begin{figure*}[t]
	\centering
	\subfigure[Original Image.]{
		\label{SPEC1}

		\includegraphics[height=2.3cm,width=2.4cm]{blood_ori.jpg}}
	\hfill
	\subfigure[Superpixel.]{
		\label{SPEC2}
	    \includegraphics[height=2.3cm,width=2.4cm]{blood_super.jpg}}
	    \hfill
	\subfigure[K-Means.]{
		\label{CELL1}

		\includegraphics[height=2.3cm,width=2.7cm]{KM1.jpg}}
	\hfill
	\subfigure[GSP.]{
		\label{CELL2}
		\includegraphics[height=2.3cm,width=2.7cm]{GSP1.jpg}}
	\hfill
	\subfigure[MLN-FZ.]{
		\label{CELL3}
		\includegraphics[height=2.3cm,width=2.7cm]{MLN_CP1.jpg}}
	\hfill
	\subfigure[MLN-SVD.]{
		\label{CELL4}
		\includegraphics[height=2.3cm,width=2.7cm]{MLN_HOSVD.jpg}}
	\vspace{-2mm}
	\caption{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).}
	\vspace{-4mm}
	\label{SPEC}
\end{figure*}
\section{Application Examples} \label{app}
We now provide some illustrative application examples within our
M-GSP framework.

\vspace{-3mm}
\subsection{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. \ref{SPEC2}. 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 \ref{basic4}. 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. (\ref{decompose1}), and
singular vector in Eq. (\ref{decomposeS}).
Step 5 determines
$K$ based on the largest arithmetic
gap in eigenvalues.
Traditional clustering methods, such as $k$-means clustering \cite{d11}, can be carried out in Step 6.

\begin{table}[t]
	\caption{Results of Image Segmentation in Image BSD300: (F) for all, and (C) for coarse}
	\vspace{-2mm}
	\scriptsize
	\begin{tabular}{|l|l|l|l|l|l|}
		\hline
		& N=100 (F)  & N=300 (F)   & N=100 (C) & N=300 (C) & N=900 (C)   \\ \hline
		GSP       & 0.1237          & 0.1149          & 0.3225                & 0.3087                & 0.3067                      \\ \hline
		K-MEANS   & 0.1293          & 0.1252          & 0.3044                & 0.3105                & 0.3124                        \\ \hline
		MLN-SVD & \textbf{ 0.1326}    & \textbf{0.1366} & \textbf{ 0.3344}          & \textbf{0.3394}     & \textbf{ 0.3335}         \\ \hline
		MLN-CP    & 0.1321          & 0.1293          & 0.3195                & 0.3256                & 0.3243                      \\ \hline
	\end{tabular}
	\label{BSD}
\end{table}

\begin{algorithm}[t]
	\begin{algorithmic}[1]
		\caption{MLN-based Unsupervised Image Segmentation}\label{basic4}
		\STATE {\bf{Input}}: RGB Image $\mathbf{I}\in\mathbb{R}^{P\times Q\times 3}$;
		\STATE 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}$;
		\STATE Construct a multilayer network $\mathbf{A}\in\mathbb{R}^{M\times N \times M \times N}$;
		\STATE Find entity-wise spectrum $\mathbf{E}=[\mathbf{e}_1,\cdots,\mathbf{e}_N]\in\mathbb{R}^{N\times N}$;
		\STATE 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}$;
		\STATE 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;
		\STATE Assign all pixels inside one superpixel to the cluster of that superpixel;
		\STATE  {\bf{Output}}: The segmented image.
	\end{algorithmic}
\end{algorithm}

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. \ref{SPEC1}.
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 \ref{CELL1}-\ref{CELL4}.
{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 \textbf{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 \cite{d13}. 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 \ref{BSD}, 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.

\subsection{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$
\begin{equation}
\min_\mathbf{a}\quad ||M(f(\mathbf{s}))_L-\mathbf{y}_L||_2^2,
\end{equation}
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.

\begin{figure}[t]
	\centering
	\includegraphics[width=2in]{class3.jpg}
		\vspace{-2mm}
	\caption{Results of Classification.}
	\label{cl}
	\vspace{-5mm}
\end{figure}

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
\begin{equation}\label{sss}
\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},
\end{equation}
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. (\ref{sss}) 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}'=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. \ref{cl},
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.

\vspace{-3mm}
\subsection{Dynamic Point Cloud Analysis}
\begin{figure}[t]
	\centering
	\includegraphics[width=3in]{Figure2.png}
	\vspace{-2mm}
	\caption{Example of Transformed Signals in a Dynamic Point Cloud.}
	\label{ex}
	\vspace{-5mm}
\end{figure}

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 \cite{CMU}. 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$ \cite{c6}.
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(-\Vert \mathbf{x}_i-\mathbf{x}_j \Vert^2_2/{\sigma^2}
)$ only if $\Vert \mathbf{x}_i-\mathbf{x}_j \Vert^2_2\leq\tau_s$.
Next, we estimate the spatial and temporal basis vectors of each shorter-term sequences by HOSVD in Eq. (\ref{decomposeS}). 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) \cite{c36}.
In Fig. \ref{ex}, 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. \ref{ex}, 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 \textbf{Appendix F}.

\section{Conclusion} \label{con}
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.
\IEEEpeerreviewmaketitle

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi