# Time-Varying Graph Mode Decomposition

Naveed ur Rehman  N. Rehman is with the Department of Electrical and Computer
Engineering, Aarhus University, 8200 Aarhus N, Denmark e-mail:
<EMAIL_ADDRESS>
###### Abstract

Time-varying graph signals are alternative representation of multivariate (or
multichannel) signals in which a single time-series is associated with each of
the nodes or vertex of a graph. Aided by the graph-theoretic tools, time-
varying graph models have the ability to capture the underlying structure of
the data associated with multiple nodes of a graph – a feat that is hard to
accomplish using standard signal processing approaches. The aim of this
contribution is to propose a method for the decomposition of time-varying
graph signals into a set of graph modes. The graph modes can be interpreted in
terms of their temporal, spectral and topological characteristics. From the
temporal (spectral) viewpoint, the graph modes represent the finite number of
oscillatory signal components (output of multiple band-pass filters whose
center frequencies and bandwidths are learned in a fully data-driven manner),
similar in properties to those obtained from the empirical mode decomposition
and related approaches. From the topological perspective, the graph modes
quantify the functional connectivity of the graph vertices at multiple scales
based on their signal content. In order to estimate the graph modes, a
variational optimization formulation is designed that includes necessary
temporal, spectral and topological requirements relevant to the graph modes.
An efficient method to solve that problem is developed which is based on the
alternating direction method of multipliers (ADMM) and the primal-dual
optimization approach. Finally, the ability of the method to enable a joint
analysis of the temporal and topological characteristics of time-varying graph
signals, at multiple frequency bands/scales, is demonstrated on a series of
synthetic and real time-varying graph data sets.

###### Index Terms:

graph signal processing, multivariate signal decomposition, network analysis,
multiscale connectivity networks

##  I Introduction

Graph signal processing (GSP) is an emerging multidisciplinary research area
that deals with the concepts, theory and methods for studying signals residing
on irregular structures or networks [1]. Signals exhibiting those
characteristics are now routinely collected from diverse sources and sensors
in modern networks e.g., biological [2], energy [3], sensor- and social-
networks [4]. Classical signal processing approaches, mainly designed for
signals residing in the euclidean spaces (time-series and images), fall short
of effectively modeling the signals defined on irregular structures and
elucidating their complex interactions.

To this end, GSP aims to extend classical signal processing tools to signals
supported on irregular (graph) structures. This has led to new concepts and
methods extending the harmonic analysis to graph signals e.g., graph Fourier
transform [5]; filter bank [6], compressed sensing and reconstruction of graph
signals [7]; estimating graph spectrum through related new concepts of
stationarity for graphs [8]; and vertex-frequency graph representations [9].

While the above concepts and methods are mainly applicable to static graph
signals i.e., those that do not evolve with time, recent efforts are also
targeted towards the analysis of graph signals that are dynamic (or time-
varying) in nature. Those are referred to as time-varying graph signals. In
essence, time-varying graph signals are intimately related to multivariate (or
multichannel) signals that comprise time-series data associated with its each
variate (or channel). Using graph-theoretic tools, such data could be
alternatively modeled as time-varying graph signals in which each time-series
is associated with a node (or a vertex) of a graph, with the topological
interactions between the graph nodes being represented through an adjacency
(or Laplacian) matrix. Indeed, the ability of graph-theoretic framework to
effectively model the nodal connections within irregularly structured data
makes it a more suitable candidate to handle multivariate signals as compared
to the standard signal processing approaches. Some notable examples of graph-
powered techniques for multivariate signals include frameworks for harmonic
time-vertex analysis [10], time-varying graph reconstruction [11], subsampling
[12] and denoising [13]. Still, graph-theoretic tools are largely an untapped
resource when it comes to information extraction from multivariate signals.

Over the last couple of decades, one of the areas in non-stationary signal
processing that has found interdisciplinary applications is the data-driven
harmonic decomposition of non-stationary signals e.g., empirical mode
decomposition (EMD) [14], variational mode decomposition (VMD) [15],
synchrosqueezed transform (SST) [16] and sliding-window singular spectrum
analysis (SSA) [17]. Originally designed for single-channel time series, the
overarching goal of signal decomposition methods is to decompose a signal into
its weakly non-stationary and oscillatory amplitude- frequency-modulated (AM-
FM) components in a fully data-driven manner i.e., without prior assumptions
on input data. Recently, extensions of these methods for multivariate signals
have also emerged, namely multivariate extensions of empirical mode
decomposition (MEMD) [18], variational mode decomposition (MVMD) [19] and
iterative filtering [20]; a recent article provides an overview of both
univariate and multivariate signal decomposition methods along with their
comparative analysis [21]. As those are derived directly from their univariate
counterparts that focused only on the temporal variations of a signal by
design, the existing multivariate extensions fail to model complex structural
(inter-channel) connections within data. As a result, none of the existing
multivariate decomposition approaches are able to provide information about
the rich functional connectivity patterns that maybe present in the data. In
this paper, this issue is addressed through the unification of the concepts
from GSP and harmonic analysis by modeling multivariate signals as time-
varying graphs.

Specifically, this contribution proposes the first-ever model for the
decomposition of time-varying graph signals into multiple graph modes using
variational optimization. The developed model includes information about the
network (or topological) interactions as well as the dynamic (time-varying)
behaviour of the data. The desired characteristics of the graph modes dictates
the design of the resulting optimization formulation - reflecting both the
temporal and spectral requirements as well as the topological constraints that
the graph modes must exhibit. For instance, the temporal (and spectral)
requirements include the separation of the input signal into a finite number
of oscillatory (and band-limited) components that are derived in a fully data-
driven manner i.e., without making any assumptions in the form of basis
functions. Those components must also sum up to obtain the original signal
(reconstruction). From the topological perspective, the following
characteristics of the graph modes are reflected in the optimization
formulation: i) the modes should include explicit information about the nodal
connections within data, via associated adjacency or Laplacian matrix; ii) the
modes should be smooth over the respective graph structure. We employ the
alternating direction method of multipliers (ADMM) and primal-dual approach to
solve the resulting optimization formulation that alternates between the
estimation of the band-limited signal components of the graph modes and the
estimation of their respective topological structures. The resulting method,
termed the time-varying graph mode decomposition (TVGMD), enables the joint
analysis of signal’s temporal dynamics and the corresponding network
structures at multiple scales.

##  II Related Work

Recently, temporal graph signal decomposition (TGSD) method has been proposed
for the decomposition of time-varying graph signals [22]. The key idea behind
this method is to encode an input signal through a combination of the fixed
graph and fixed time dictionaries, resulting in a compact representation of
multivariate data. In essence, the TGSD method cleverly extends the recent
graph dictionary learning approaches to account for the decomposition of time-
varying signals [23]. The proposed method (TVGMD), however, is fundamentally
different and superior to the TGSD method in the following ways: i) TVGMD is
fully data-driven, learning the inherent time-varying signal components as
well as the associated connectivity graphs directly from the data, whereas the
TGSD method uses over-complete dictionaries (or basis functions), making it a
projection-based method; ii) TVGMD provides explicit information regarding the
data connectivity of the graph nodes at multiple scales, which TGSD does not
accomplish.

A different class of multivariate signal decomposition methods are inspired
from univariate data-driven signal decomposition approaches [18, 19, 20]. To
reiterate, this family of methods focus mainly on the temporal features of a
signal, ignoring the rich topological structures that may exist between
multiple data sources (channels). For instance, the only topological
constraint included within the MVMD method is the enforcement of the same
center frequency across all data channels, in a single output mode (or
component) [19]. Since the center frequency is a global signal property, the
MVMD (and other decomposition methods) fail to model, enforce or capture local
signal characteristics e.g., signal smoothness across connected data channels
at multiple scales. Further, no direct mechanism exists to provide information
regarding the functional connectivity of data by using these approaches.

By modeling the multivariate signals as time-varying graph signals and
designing a variational optimization formulation that includes both the
temporal and topological constraints for graph modes, the proposed TVGM method
accomplishes a fully data-driven decomposition of time-varying graph signals.
The method concurrently obtains multiple graph signal modes and the associated
connectivity graphs which both the TGSD and multivariate decomposition methods
fail to achieve.

The rest of the paper is organized as follows: Section III reviews the
preliminary concepts of graph signal processing (GSP), introduces the problem
statement for TVGMD and specifies the associated optimization formulation.
Section IV elucidates the steps to solve the TVGMD optimization problem.
Section V demonstrates the utility of the TVGMD method through detailed
experiments on a range of synthetic and real-life time-varying graph signals.
Specifically, the real time-varying graph signals used in our experiments
include electroencephalogram (EEG), electricity consumption signals and plant-
wide oscillation data from an industrial setting. The paper concludes with the
discussion and conclusion sections.

##  III Time-varying Graph Mode Decomposition: Setting up the Optimization
formulation

###  III-A Preliminaries on Graph Signals

We denote a weighted and undirected graph by $R=(\mathcal{V},\mathcal{E},W)$,
where $\mathcal{V}$ represents the finite set of $N$ nodes or vertices and
$\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$ denotes the finite set of
edges of the graph. $W\in\mathbb{R}^{N\times N}$ is the weighted and symmetric
adjacency matrix of the graph. The entry $W_{mn}$ of the weighted adjacency
matrix $W$ quantifies the weight of the edge between the $m$-th and the $n$-th
vertices. These weights are generally non-negative and greater the weight,
greater the connectivity (or similarity) between the corresponding vertices.
$W_{mn}=0$ means that the vertices $m$ and $n$ are not connected.

Another matrix that is widely used to represent the graph connectivity is the
graph Laplacian matrix, denoted by $L$. It can be defined in terms of $W$ as
$L=D-W$, where $D=$ diag$(d_{1},d_{2},\ldots,d_{N})$ denotes the degree matrix
that is formed from the vertex degrees. Particularly, $d_{n}$, which is the
$n$-th diagonal entry of $D$, is computed by taking the sum of weights of all
the edges that are connected to the $n$-th vertex.

We define a graph signal on $R$ as a function $u:\mathcal{V}\to\mathbb{R}$
which assigns a real (signal) value to each vertex. This way, a graph signal
can be represented as a $N$-dimensional vector $\mathbf{u}\in\mathbb{R}^{N}$
whose $n$-th entry, $u_{n}$, denotes the graph signal measurement (or value)
at the $n$-th vertex.

In a similar vein, a time-varying graph signal on $R$ can be defined as a
function $\mathbf{u}:\mathcal{V}\to\mathbb{R}^{T}$ that assigns a (time-
series) vector of length $T$, to each vertex. With a time-series assigned to
each of the $N$ nodes of the graph, the resulting time-varying graph signal
can be represented as a matrix $U\in\mathbb{R}^{N\times T}$. The $t$-th column
of $U$, denoted by $\mathbf{u}_{t}$, represents the graph signal corresponding
to the $t$-th time instant. Further, the $n$-th row of $U$, denoted by
$\tilde{\mathbf{u}}_{n}$, represents the time-series corresponding to the
$n$-th vertex of the graph $R$. In terms of notation for multivariate signal,
$\tilde{\mathbf{u}}_{n}$ represents the $n$-th channel of the multivariate
signal $U$. Note that in the limiting case of $T\to\infty$ and infinite
sampling frequency, a continuous and unlimited-duration time-varying graph
signal $\tilde{u}_{n}(t)$ will replace $\tilde{\mathbf{u}}_{n}$ \-
representing a continuous-time function rather than a vector, corresponding to
the $n$-th vertex.

The smoothness of a time-varying graph signal $U$ on $R$ can be measured and
quantified in terms of the following function of the graph Laplacian $L$

$\displaystyle Tr\big{(}U^{T}LU\big{)}$
$\displaystyle=\frac{1}{2}\sum_{m}\sum_{n}W_{mn}||\tilde{\mathbf{u}}_{m}-\tilde{\mathbf{u}}_{n}||^{2}$
(1) $\displaystyle=\frac{1}{2}||W\circ Z||_{1,1},$

where $W_{mn}\in\mathbb{R}_{+}$ denotes the weight of the edge between the
nodes $m$ and $n$. the resulting matrix $W\in\mathbb{R}^{N\times N}$ is the
weighted graph adjacency matrix that is related to the graph Laplacian matrix
via $L=D-W$. Further,
$Z_{mn}=||\tilde{\mathbf{u}}_{m}-\tilde{\mathbf{u}}_{n}||^{2}$ are the
elements of the pairwise distance matrix $Z\in\mathbb{R}_{+}^{N\times N}$,
$\circ$ is the Hadamard product, $Tr(.)$ denotes the trace of a matrix and
$||Z||_{1,1}$ is the element-wise norm-1 of $Z$. It can be verified that the
smaller the value of the above measure (1), the greater the smoothness of the
graph signal over $R$ i.e., strongly connected vertices (with higher
corresponding edge weights) will have similar signal values.

###  III-B Signal Decomposition into Band-limited Modes

Signal decomposition (SD) refers to the process of decomposing a uni- or
multi-variate signal into its constituent finite number of band-limited modes
or components. Let $\tilde{x}(t)$ and $\tilde{u}_{k}(t)$ denote the univariate
input signal and its extracted $k$-th mode respectively, then SD results in
the following

$\tilde{x}(t)=\sum_{k=1}^{K}\tilde{u}^{(k)}(t).$ (2)

Several data-driven approaches to SD have emerged in recent years that
accomplish this decomposition without making assumptions on data e.g., without
basis functions or dictionaries. One of those methods is variational mode
decomposition (VMD) [15] that uses a variational optimization program to
obtain the desired modes. The cost function and the constraints of the
optimization program respectively are the sum of the bandwidths of all modes
${\\{}\tilde{u}^{(k)}(t){\\}}_{k=1}^{K}$ and the reconstruction of the signal
from its modes. The optimization program for VMD is as follows

$\displaystyle\underset{\\{\tilde{u}^{(k)}\\},\\{\omega^{(k)}\\}}{\text{\normalsize{minimize}}}$
$\displaystyle\Bigg{\\{}{\sum_{k}}\Bigg{\|}{\partial_{t}\Big{[}\tilde{u}_{+}^{(k)}(t)e^{-j\omega^{(k)}t}\Big{]}\Bigg{\|}^{2}_{2}}\Bigg{\\}}$
(3) subject to $\displaystyle\sum_{k}\tilde{u}^{(k)}(t)=\tilde{x}(t),$

where $u_{+}^{(k)}(t)$ denotes an analytic signal corresponding to
$u^{(k)}(t)$; the symbol $\partial_{t}$ represents partial derivative
operation with respect to time; $\\{u^{(k)}\\}$ and $\\{\omega^{(k)}\\}$
respectively denote sets of all $K$ number of modes and their center
frequencies. The squared $l^{2}$-norm term within the parenthesis in the cost
function (3) is an estimator of the sum of the bandwidths of all the signal
modes.

An extension of VMD to multivariate signals was developed in [19] that
utilized a model for multivariate modulated oscillations to design the
following optimization formulation for the decomposition of multivariate
signals.

$\displaystyle\underset{\\{u_{n}^{(k)}\\},\\{\omega^{(k)}\\}}{\text{minimize}}$
$\displaystyle\Bigg{\\{}{\sum_{k}\sum_{n}}\Bigg{\|}\partial_{t}\Big{[}\tilde{u}_{n+}^{(k)}(t)e^{-j\omega^{(k)}t}\Big{]}\Bigg{\|}^{2}_{2}\Bigg{\\}}$
(4) subject to
$\displaystyle\sum_{k}\tilde{u}_{n}^{(k)}(t)=\tilde{x}_{n}(t),\text{
}n=1,2,\ldots,N.$

Note that the main focus of MVMD is on extracting the temporal features of a
signal i.e., band-limited oscillations. The only topological constraint that
it enforces on a signal mode is to have the same fixed center frequency
$\omega^{(k)}$ across all its channels. Clearly, that is inadequate to model
and extract rich topological structures that time-varying graph signals
typically possess.

###  III-C Problem Statement

Consider a time-varying graph signal $\mathbf{\tilde{x}(t)}$ that comprises a
set of time-series signals:
$\mathbf{\tilde{x}(t)}=[\tilde{x}_{1}(t),\ldots,\tilde{x}_{N}(t)]$, with
$\tilde{x}_{n}(t)$ representing the time-series corresponding to the $n$-th
vertex. In the discrete domain, $\mathbf{\tilde{x}(t)}$ could also be
represented by a matrix $\mathbf{X}\in\mathbb{R}^{N\times T}$ where $T$
denotes the number of time samples of the time-varying graph signal.

Let the graph signal $\mathbf{\tilde{x}(t)}$ be decomposed into a set of band-
limited oscillatory components as follows

$\mathbf{\tilde{x}}(t)=\sum_{k=1}^{K}\mathbf{g}^{(k)}(t),$ (5)

where
$\mathbf{g}^{(k)}(t)=[g_{1}^{(k)}(t),g_{2}^{(k)}(t),\ldots,g_{N}^{(k)}(t)]$.
The set of graph modes of $\mathbf{\tilde{x}(t)}$ comprise its band-limited
constituent components along with their weighted adjacency or connectivity
matrices. For instance, the $k$-th graph mode is represented by
$G^{(k)}=[\mathbf{g}^{(k)}(t),W^{(k)}]$, with the collection of all the graph
modes of $\mathbf{\tilde{x}(t)}$ be denoted by
$\mathbf{G}=[G^{(1)},G^{(2)},\ldots,G^{(K)}]$. Concretely, given the number of
modes $K$ as defined by the user, the problem here is to decompose a time-
varying graph signal $\mathbf{\tilde{x}(t)}$ into its constituent graph modes

$\mathbf{\tilde{x}(t)}\to[G^{(1)},G^{(2)},\ldots,G^{(K)}].$ (6)

###  III-D Optimization Formulation for Time-Varying Graph Mode Decomposition
(TVGMD)

To obtain the desired graph modes $\mathbf{G}$ from input time-varying graph
signal $\mathbf{\tilde{x}(t)}$, we design a variational optimization program
that includes both the temporal and the topological constraints for the
extraction of physically meaningful modes. Specifically, the temporal
requirements include the extraction of AM-FM oscillatory components,
$\mathbf{g}^{(k)}(t)=[g_{1}^{(k)}(t),g_{2}^{(k)}(t),\ldots,g_{N}^{(k)}(t)]$,
which are band-limited in nature, with the center frequency and the bandwidth
of each mode being estimated via a data-driven approach. Another temporal
requirement is the exact reconstruction of the original graph signal from its
extracted modes.

Next, topological or structural constraints related to the time-varying graph
modes have been included in the optimization formulation. Particularly, in the
proposed optimization program, we embed the requirement of the smoothness of
each extracted mode $\mathbf{g}^{(k)}(t)$ with respect to its graph structure
$W^{(k)}$. To this end, the (non-)smoothness measure (1) corresponding to the
$K$ graph modes $\\{G^{(k)}\\}_{k=1}^{K}$ is included in the optimization
program. Finally, the optimization program includes requirements for the
validity of the graph adjacency matrices. The resulting optimization
formulation for time-varying graph mode decomposition (TVGMD) method is given
as follows

$\displaystyle\begin{split}\underset{\\{\tilde{g}_{n}^{(k)}\\},\\{\omega^{(k)}\\},\\{W^{(k)}\\}}{\text{minimize}}\quad&\Bigg{\\{}\alpha{\sum_{k=1}^{K}\sum_{n=1}^{N}}\Big{\|}\partial_{t}\Big{[}\tilde{g}_{n+}^{(k)}(t)e^{-j\omega^{(k)}t}\Big{]}\Big{\|}^{2}_{2}\Bigg{\\}}\\\
&\quad+\beta{\sum_{k=1}^{K}}||W^{(k)}\circ Z^{(k)}||_{1,1}\\\
&\quad+\frac{\gamma}{2}{\sum_{k=1}^{K}}\lVert
W^{(k)}\rVert^{2}_{F}-{\sum_{k=1}^{K}}\mathbf{1}^{T}\log(W^{(k)}\mathbf{1})\end{split}$
(7) subject to
$\displaystyle\sum_{k}\tilde{g}_{n}^{(k)}(t)=\tilde{x}_{n}(t),\quad\forall n$
(8) $\displaystyle W^{(k)}\geq\mathbf{0},\quad W^{(k)}=W^{(k)^{T}},$
$\displaystyle\text{diag}(W^{(k)})=0,\quad\forall k$

where $\tilde{g}_{n+}^{(k)}(t)$ denotes an analytical signal corresponding to
$\tilde{g}_{n}^{(k)}(t)$. The output of the optimization program includes the
graph modes $\\{G^{(k)}\\}_{k=1}^{K}$ that comprise the time-varying
oscillatory graph components, $\\{\mathbf{g}^{(k)}\\}_{k=1}^{K}$, along with
their respective connectivity structures $\\{W^{(k)}\\}_{k=1}^{K}$. The center
frequency $\omega^{(k)}$ of each of the extracted graph modes is also
estimated.

Some comments on the TVGMD optimization formulation are in order. The first
term of the cost function (7) represents the sum of the the bandwidths of all
the extracted graph modes that ensures the extraction of band-limited
oscillatory components in the temporal domain [19]. The second term in the
cost function imposes the smoothness of the extracted graph modes over their
respective graph structures $W^{(k)}$. This way, valid graph structures
$W^{(k)}$ are obtained i.e., for a pair of nodes with similar signal values,
the corresponding adjacency matrix element has a higher weight. Further, the
term involving the Frobenius norm of $W^{(k)}$ prohibits the formation of very
large edge weights in $W^{(k)}$ while not penalizing the smaller weights. The
last term can be interpreted as a logarithmic barrier on the node degree
vector, enforcing the positive degrees on the graph nodes while not preventing
the individual connections from becoming zero; this results in a better
connected graph [24].

As for the constraints of the optimization problem, the first constraint
guarantees that the obtained graph modes fully reconstruct the original graph
signal. The remaining constraints relate to the validity of the adjacency
matrices at multiple scales, including their symmetricity, non-negative values
and a zero vector at the diagonal.

##  IV Optimization Solution

The first step towards the solution of the TVGMD optimization problem is to
convert the original constrained optimization formulation to the unconstrained
formulation. Then, a combination of the alternating direction method of
multipliers (ADMM) and the primal-dual approach is used to solve the problem.

Note that in the original optimization problem there are constraints related
to the symmetricity of $W^{(k)}$ along with a zero vector as its diagonal.
Those constraints are cumbersome while searching for a convenient solution to
the problem. To address this problem, it was proposed in [24] to use a vector
representation of $W^{(k)}$ within the space
$\mathcal{W}_{v}=w\in\mathbb{R}_{+}^{N(N-1)/2}$. Following that approach, let
$w^{(k)}$ and $z^{(k)}$ denote the vectors comprising all entries above the
main diagonal of $W^{(k)}$ and $Z^{(k)}$ respectively. Further, $Q^{(k)}$
denotes a binary matrix that fulfills the condition
$Q^{(k)}w^{(k)}=W^{(k)}\mathbf{1}$. This way,
$Q^{(k)}w^{(k)}\in\mathbb{R}^{N}$ has a convenient interpretation as a vector
of node degrees for $\mathbf{g^{(k)}}(t)$. Then, the original optimization
problem can be modified as follows

$\displaystyle\begin{split}\underset{\\{\tilde{g}_{n}^{(k)}\\},\\{\omega^{(k)}\\},\\{w^{(k)}\\}}{\text{minimize}}\quad&\Bigg{\\{}\alpha{\sum_{k=1}^{K}\sum_{n=1}^{N}}\Big{\|}\partial_{t}\Big{[}\tilde{g}_{n+}^{(k)}(t)e^{-j\omega^{(k)}t}\Big{]}\Big{\|}^{2}_{2}\Big{\\}}\\\
&\quad+\beta{\sum_{k=1}^{K}}2w^{(k)T}z^{(k)}\\\
&\quad+\gamma{\sum_{k=1}^{K}}\lVert
w^{(k)}\rVert^{2}_{2}-{\sum_{k=1}^{K}}\mathbf{1}^{T}\log(Q^{(k)}w^{(k)})\end{split}$
(9) subject to
$\displaystyle\sum_{k}\tilde{g}_{n}^{(k)}(t)=\tilde{x}_{n}(t),\quad\forall n$
(10) $\displaystyle w^{(k)}\geq\mathbf{0}.\quad\forall k$

By introducing the sets of dual variables $\lambda_{n}(t)$, corresponding to
each of the equality constraints in (10), the augmented Lagrangian function of
the above optimization formulation can be written as

$\displaystyle\mathcal{L}(\\{\tilde{g}_{n}^{(k)}\\},\\{\omega^{(k)}\\},\\{w^{(k)}\\})=h_{1}(\\{\tilde{g}_{n}^{(k)}\\},\\{\omega^{(k)}\\})+h_{2}(\\{w^{(k)}\\})$
(11)

where

$h_{1}(\\{\tilde{g}_{n}^{(k)}\\},\\{\omega^{(k)}\\})=\alpha{\sum_{k=1}^{K}\sum_{n=1}^{N}}\Big{\|}\partial_{t}\Big{[}\tilde{g}_{n+}^{(k)}(t)e^{-j\omega^{(k)}t}\Big{]}\Big{\|}^{2}_{2}+{\sum_{n}}\Big{\|}\tilde{x}_{n}(t)-\sum_{k}\tilde{g}_{n}^{(k)}(t)\Big{\|}_{2}^{2}+{\sum_{n}}\Big{\langle}\lambda_{n}(t),\tilde{x}_{n}(t)-\sum_{k}\tilde{g}_{n}^{(k)}(t)\Big{\rangle},$
(12)
$h_{2}(\\{w^{(k)}\\})=\beta{\sum_{k=1}^{K}}2w^{(k)T}z^{(k)}+\gamma{\sum_{k=1}^{K}}\lVert
w^{(k)}\rVert^{2}_{2}-{\sum_{k=1}^{K}}\mathbf{1}^{T}\log(Q^{(k)}w^{(k)})+{\sum_{k=1}^{K}}\mathbb{I}_{\\{w^{(k)}\geq
0\\}}$ (13)

with $\mathbb{I}_{\\{w\geq 0\\}}$ representing the indicator function that is
equal to zero for $w\geq 0$ but is otherwise infinite. The term involving the
indicator function has been added to cater for the inequality constraints in
(10). Note the difference between the symbols $\omega^{(k)}$ and $w^{(k)}$
that respectively denote the center frequency and the weight vector
(corresponding to $W^{(k)}$) of the $k$-th graph mode.

Algorithm 1 ADMM Optimization Strategy for TVGMD

Initialize:

$\\{\tilde{g}_{n,1}^{(k)}\\}\leftarrow 0\mbox{,
}\\{\omega^{(k)}_{1}\\}\leftarrow 0\mbox{, }\\{w^{(k)}_{1}\\}\leftarrow
0\mbox{, }\lambda_{n,1}\mbox{, }i\leftarrow 0\mbox{, }\epsilon\leftarrow
10^{-7}$

repeat

$i\leftarrow i+1$

for $k=1:K$ do

for $n=1:N$ do Update graph modes in time $\tilde{g}_{n}^{(k)}$:

$\tilde{g}_{n,i+1}^{(k)}/\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k)}\leftarrow\underset{\tilde{g}_{n}^{(k)}}{\text{arg min
}}\mathcal{L}\left({\\{}\tilde{g}_{n,i}^{(k^{\prime})}{\\}},{\\{}\omega_{i}^{(k^{\prime})}{\\}},{\\{}w_{i}^{(k^{\prime})}{\\}},\lambda_{n,i}\right)$
(14)

end for

end for

for $k=1:K$ do Update center frequency $\omega^{(k)}$:

$\omega_{i+1}^{(k)}\leftarrow\underset{\omega^{(k)}}{\text{arg min
}}\mathcal{L}\left({\\{\tilde{g}_{n,i+1}^{(k^{\prime})}\\},\\{\omega_{i}^{(k^{\prime})}\\},\\{w_{i}^{(k^{\prime})}\\},\lambda_{n,i}}\right)$
(15)

end for

for $k=1:K$ do Update modes along graph geodesics:

$U_{i+1}^{(k)}\leftarrow\underset{\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern 1.5mu^{(k)}}{\text{arg min
}}\mathcal{L}\left({\\{}\mkern 1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k^{\prime})}{\\}},{\\{}\omega_{i+1}^{(k^{\prime})}{\\}},{\\{}w_{i}^{(k^{\prime})}{\\}},\lambda_{n,i}\right)$
(16)

end for

for $k=1:K$ do Update graph adjacency matrix $w^{(k)}$:

$w_{i+1}^{(k)}\leftarrow\underset{w^{(k)}}{\text{arg min
}}\mathcal{L}\left({\\{\tilde{g}_{n,i+1}^{(k^{\prime})}\\},\\{\omega_{i+1}^{(k^{\prime})}\\},\\{w_{i}^{(k^{\prime})}\\},\lambda_{n,i}}\right)$
(17)

end for

for $n=1:N$ do Update $\lambda_{n}$:

$\lambda_{n,i+1}=\lambda_{n,i}+\tau\Big{(}\tilde{x}_{n}-\sum_{k}\tilde{g}_{n,i+1}^{(k)}\Big{)}$
(18)

end for

until Convergence:
$\sum_{k}\sum_{n}\frac{\|\tilde{g}_{n,i+1}^{(k)}-\tilde{g}_{n,i}^{(k)}\|_{2}^{2}}{\|\tilde{g}_{n,i}^{(k)}\|_{2}^{2}}<\epsilon$

The unconstrained optimization problem (11) is solved using the alternating
direction method of multipliers (ADMM) [25] and the primal-dual approach [26].
The ADMM operates by converting a complex optimization problem into a series
of simpler sub-optimization problems. Our approach towards the solution of
(11) is outlined in the Algorithm 1. By using the ADMM approach, the TVGMD
unconstrained optimization problem (11) is transformed into 4 simple
optimization problems, relating to the: i) update of the graph modes
$\tilde{g}_{n}^{(k)}$ in the time domain (14); ii) update of the center
frequency $\omega^{(k)}$ for each graph mode (15); iii) update of the graph
modes along the graph geodesics or vertices $U^{(k)}$ (16); and iv) update of
the adjacency or Laplacian matrix $W^{(k)}$ for each graph mode (17). Finally,
the dual variables $\lambda_{n}$ are updated according to (18). In the
following, we obtain the solutions to the above optimization sub-problems.

###  IV-A Graph mode update in the time domain:

To solve the minimization problem related to the graph mode update in the time
domain (14), we write the associated optimization formulation below

$\displaystyle\tilde{g}_{n,i+1}^{(k)}=\underset{\tilde{g}_{n}^{(k)}}{\text{arg
min}}\Bigg{\\{}\alpha\Big{\|}\partial_{t}\Big{[}\tilde{g}_{n+}^{(k)}(t)e^{-j\omega^{(k)}t}\Big{]}\Big{\|}^{2}_{2}+$
(19)
$\displaystyle\Big{\|}\tilde{x}_{n}(t)-{\sum_{k^{\prime}}}\tilde{g}_{n}^{(k^{\prime})}(t)+\frac{\lambda_{n}(t)}{2}\Big{\|}^{2}_{2}\Bigg{\\}}.$

The above problem was encountered within the MVMD optimization formulation
[19] and was solved in the Fourier (or frequency) domain. Using that result,
the following graph mode update relation (in the time domain) is obtained

$\displaystyle\hat{g}_{n,i+1}^{(k)}(\omega)=\frac{\hat{x}_{n}(\omega)-\sum_{k^{\prime}\neq
k}\hat{g}_{n}^{(k^{\prime})}(\omega)+\frac{\hat{\lambda}_{n}(\omega)}{2}}{1+2\alpha(\omega-\omega^{(k)})^{2}},$
(20)

where $\hat{g}_{n}^{(k)}(\omega)$ denotes the Fourier transform of the $k$-th
graph mode, $\tilde{g}_{n}^{(k)}(t)$.

###  IV-B Center frequency update:

The center frequency of each graph mode is updated by solving the optimization
problem in (15). Here, only the first term of (12) depends on $\omega_{k}$,
simplifying (15) to

$\omega_{n,i+1}^{(k)}=\underset{\omega^{(k)}}{\text{arg
min}}\Big{\\{}{\sum_{n}}\Big{\|}\partial_{t}\Big{[}\tilde{g}_{n+}^{(k)}(t)e^{-j\omega^{(k)}t}\Big{]}\Big{\|}^{2}_{2}\Bigg{\\}}.$
(21)

Like the mode update problem, the above optimization problem was also solved
in the frequency domain within the MVMD optimization solution; see (30)-(32)
in [19]. The resulting update relation for the center frequency $\omega^{(k)}$
of the $k$-th graph mode is given by

$\omega_{i+1}^{(k)}=\frac{{\sum_{n}}\int^{\infty}_{0}\omega\times|\hat{g}_{n}^{(k)}(\omega)|^{2}d\omega}{{\sum_{n}}\int^{\infty}_{0}|\hat{g}_{n}^{(k)}(\omega)|^{2}d\omega}.$
(22)

###  IV-C Graph mode update along the graph geodesics:

The next step involves the graph mode update along the graph geodesics (or
vertices). This is fundamentally different from the mode update in the time
domain which was conveniently performed through the Fourier transformation via
(20). In the time domain, the key requirement was the small bandwidth of each
decomposed graph mode which was imposed through the first term of (12).

On the contrary, the crucial constraint for the graph modes update along the
graph geodesics is that the obtained modes should be consistent with the
associated connectivity matrices (e.g., adjacency or the Laplacian matrix).
Simply put, the graph modes (along the graph geodesics) must be updated in
such a way that they have similar values on the pair of nodes having a strong
connection. This key requirement is implicitly encoded in the second terms of
the TVGMD cost functions in (7) and (9), and the first term of (13) that is
part of the augmented Lagrangian function (11). In those terms, $Z^{(k)}$ and
its corresponding vector $z^{(k)}$ are directly related to the graph modes
$\tilde{g}^{(k)}(t)$ via (1). That is, $Z^{(k)}$ is the pairwise distance
matrix of the graph modes at different vertices.

Based on those observations, the corresponding optimization problem for the
mode update along graph geodesics can be written as

$\tilde{u}_{n,i+1}^{(k)}=\underset{\tilde{g}_{n}^{(k)}}{\text{arg
min}}\Bigg{\\{}{\sum_{n}}\Big{\|}\tilde{x}_{n}(t)-\sum_{k^{\prime}}\tilde{g}_{n}^{(k^{\prime})}(t)\Big{\|}_{2}^{2}+2\beta
w^{(k)T}z^{(k)}\Bigg{\\}}.$ (23)

If we let $\mkern 1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k)}\\}_{k=1}^{K}$ denote the matrix representation of the graph
modes obtained by solving (14), and further by using (1), the above
optimization problem can be expressed as

$U_{i+1}^{(k)}=\underset{\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern 1.5mu^{(k)}}{\text{arg
min}}\Bigg{\\{}\Big{\|}X-{\sum_{k^{\prime}}}\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k^{\prime})}\Big{\|}_{2}^{2}+2\beta Tr\left[\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k)T}L_{i}^{(k)}\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k)}\right]\Bigg{\\}}.$ (24)

By denoting $F^{(k)}=X-\sum_{k^{\prime}\neq k}\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern 1.5mu_{i+1}^{(k^{\prime})}$,
the first term on the right hand side of (24) can be written as
$\Big{\|}F^{(k)}-\mkern 1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k)}\Big{\|}_{2}^{2}=Tr\left[(F^{(k)}-\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k)})(F^{(k)}-\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k)})^{T}\right]$, leading to the following optimization problem

$U_{i+1}^{(k)}=\underset{\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern 1.5mu^{(k)}}{\text{arg
min}}\Big{\\{}Tr\left[\Big{(}F^{(k)}-\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k)}\Big{)}\Big{(}F^{(k)}-\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k)}\Big{)}^{T}\right]+2\beta Tr\left[\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k)T}L_{i}^{(k)}\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu_{i+1}^{(k)}\right]\Big{\\}}.$ (25)

To minimize the above function that is clearly convex, we set its derivative
to zero resulting in the following update for the graph modes along the graph
geodesics

$U_{i+1}^{(k)}=\left(I+\beta L_{i}^{(k)}\right)^{-1}F^{(k)},$ (26)

where $I$ denotes the identity matrix.

###  IV-D Graph adjacency (connectivity) matrix update:

The connectivity matrices of the graph modes are updated and obtained by
solving the optimization (sub-)problem given by (17). For our unconstrained
Lagrangian function (11), the optimization problem (17) becomes

$w_{i+1}^{(k)}=\underset{w^{(k)}}{\text{arg min }}2\beta
w^{(k)T}z^{(k)}+\gamma\lVert
w^{(k)}\rVert^{2}_{2}-\mathbf{1}^{T}\log(Q^{(k)}w^{(k)})+\mathbb{I}_{\\{w^{(k)}\geq
0\\}}.$ (27)

This problem is similar in form to the popular graph learning model proposed
by Kalofolias in [24]. The only difference between the two models is the
appearance of the parameter $\beta$ in the first term of (27), which was
effectively equal to 1 in Kalofolias’ model. That problem is solved by using
the primal-dual optimization technique, illustrated in Algorithm 6 in [26]. We
use the same technique to solve (27). Particularly, to make our problem
amenable to be solved by using the primal-dual optimization technique, we cast
our optimization problem as a sum of the following three functions

$w_{i+1}^{(k)}=\underset{w^{(k)}}{\text{arg min
}}b_{1}(w)+b_{2}(Kw)+b_{3}(w),$ (28)

where

$\displaystyle b_{1}(w^{(k)})$ $\displaystyle=\mathbb{I}_{\\{w^{(k)}\geq
0\\}}+2\beta w^{(k)T}z^{(k)},$ $\displaystyle b_{2}(Q^{(k)}w^{(k)})$
$\displaystyle=-\mathbf{1}^{T}\log(Q^{(k)}w^{(k)}),$ $\displaystyle
b_{3}(w^{(k)})$ $\displaystyle=\gamma\lVert w^{(k)}\rVert^{2}.$

Using the above form of the optimization problem, the primal-dual algorithm
used for the solving that problem is given in the Appendix 1.

With the knowledge of the solutions to all the sub-optimization problems
(14)-(17) within the TVGMD algorithm, the steps to apply the TVGMD algorithm
are listed in Algorithm 2. In the algorithm, note that the TVGMD updates the
graph modes in two steps: the first update (29) relates to the requirements of
the modes in the time domain, e.g., mainly oscillatory modes having small
bandwidths, and is conducted in the Fourier domain. The second update (31)
concerns the mode constraints along the graph geodesics and is performed in
the spatial domain. Between the two updates, a conversion step (30) is needed
that transforms the modes from the frequency to the time domain in each
iteration.

Algorithm 2 Time-varying Graph Mode Decomposition

Input and user-defined parameters: $K,\alpha,\beta,\gamma\mbox{,}$

$X=\\{x_{n}(t)\\}_{n=1}^{N}\xlongleftrightarrow{FT}\\{\hat{x}_{n}(\omega)\\}_{n=1}^{N}$

Initialize:

$\\{\hat{g}_{n,1}^{(k)}\\}\leftarrow 0\mbox{,
}\\{\omega^{(k)}_{1}\\}\leftarrow 0\mbox{, }\\{W^{(k)}_{1}\\}\leftarrow
0\mbox{, }\lambda_{n,1}\mbox{, }i\leftarrow 0\mbox{, }\epsilon\leftarrow
10^{-7}$

repeat

$i\leftarrow i+1$

for $k=1:K$ do

for $n=1:N$ do Updating graph mode $\hat{g}_{n}^{(k)}$:

$\hat{g}_{n,i+1}^{(k)}(\omega)\leftarrow\frac{\hat{x}_{n}(\omega)-\sum_{k^{\prime}\neq
k}\hat{g}_{n}^{(k^{\prime})}(\omega)+\frac{\hat{\lambda}_{n}(\omega)}{2}}{1+2\alpha(\omega-\omega^{(k)})^{2}}$
(29)

end for

end for

for $k=1:K$ do Updating center frequency $\omega_{k}$:

$\omega_{n,i+1}^{(k)}\leftarrow\frac{{\sum_{n}}\int^{\infty}_{0}\omega|\hat{g}_{n,i}^{(k)}(\omega)|^{2}d\omega}{{\sum_{n}}\int^{\infty}_{0}|\hat{g}_{n,i}^{(k)}(\omega)|^{2}d\omega}$

end for

$\Big{\\{}\mkern 1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern
1.5mu^{(k)}\Big{\\}}_{k=1}^{K}\xlongleftarrow{IFT}\Big{\\{}\hat{g}_{n,i+1}^{(k)}(\omega)\Big{\\}}_{k=1}^{K}$
(30)

for $k=1:K$ do Updating modes along geodesics:

$\displaystyle L_{i}^{(k)}$ $\displaystyle\leftarrow$ $\displaystyle
D_{i}^{(k)}-W_{i}^{(k)}$ $\displaystyle U_{i+1}^{(k)}$
$\displaystyle\leftarrow$ $\displaystyle\left(I+\beta
L_{i}^{(k)}\right)^{-1}\Big{\\{}X-\sum_{k^{\prime}\neq k}\mkern
1.5mu\overline{\mkern-1.5muU\mkern-1.5mu}\mkern 1.5mu^{(k^{\prime})}\Big{\\}}$
(31)

end for

for $k=1:K$ do Updating adjacency matrix:

$\displaystyle W_{i}^{(k)}$ $\displaystyle\leftarrow$ $\displaystyle
L_{i}^{(k)}$ $\displaystyle W_{i+1}^{(k)}$ $\displaystyle\leftarrow$ Apply
Algorithm 3 in the Appendix

end for

$\Big{\\{}\hat{g}_{n,i+1}^{(k)}(\omega)\Big{\\}}_{k=1}^{K}\xlongleftarrow{FT}\Big{\\{}U_{i+1}^{(k)}\Big{\\}}_{k=1}^{K}$

for $n=1:N$ do Update $\hat{\lambda}_{c}$:

$\hat{\lambda}_{n,i+1}(\omega)=\hat{\lambda}_{n,i}(\omega)+\tau\Big{(}\hat{x}_{n}(\omega)-\sum_{k}\hat{g}_{n,i+1}^{(k)}(\omega)\Big{)}$

end for

until Convergence:
$\sum_{k}\sum_{n}\frac{\|\hat{g}_{n,i+1}^{(k)}-\hat{g}_{n,i}^{(k)}\|_{2}^{2}}{\|\hat{g}_{n,i+1}^{(k)}\|_{2}^{2}}<\epsilon$

Note that the TVGMD optimization scheme is not strictly ADMM owing to the: i)
non-convexity of the original optimization cost function (9); ii) deviation of
our optimization formulation from the standard ADMM formulation that requires
multiple set of variables having separable objective functions. Therefore, the
algorithm is not guaranteed to converge to the global minimum. That said, we
take inspiration from the fact that the ADMM approach has been successfully
employed in several popular heuristics algorithms in signal processing that
use non-convex optimization, including those based on signal decomposition
[15] and graph learning [23].

Finally, some comments on the choice of the parameters of the TVGMD algorithm
are in order. The $\alpha$ parameter determines the bandwidth of the
decomposed graph modes; the smaller the value of $\alpha$, the higher the
bandwidth of the extracted modes. For closely spaced graph modes in the
frequency domain, choosing a higher value of $\alpha$ is recommended. The
$\beta$ parameter dictates the smoothness of the graph modes
$\mathbf{g}^{(k)}(t)$ with respect to their adjacency matrices $W^{(k)}$ – an
important parameter to ensure that the adjacency (connectivity) matrices at
each scale are valid. As for the $\gamma$ parameter, it controls the sparsity
of the obtained adjacency matrices: $\gamma=0$ leads to very sparse
connectivity graphs whereas higher values tend to produce dense edge patterns
in the graphs. The parameter $\tau$ offers a trade off between the exact
reconstruction of the input graph signal (via the summation of its graph
modes) and denoising. In the presence of noise in the input signal, choosing a
very small value of $\tau$ (or even $\tau=0$) may be appropriate since exactly
reconstructing a noisy input signal will not be desirable. In section VI,
further suggestions regarding the choice of parameters in the TVGMD method are
given.

##  V Experiments, Results and Applications

This section presents the results of experiments conducted on synthetic and
real time-varying graph signals to demonstrate the power and the performance
of the TVGMD algorithm. To our knowledge, the TVGMD algorithm is the first to
concurrently obtain the constituent oscillatory components of a time-varying
graph signal and their corresponding connectivity structures in the spatial
domain, via a joint optimization scheme. To this end, the designed experiments
and the chosen applications will highlight the ability of the TVGMD algorithm
to extract both the time-domain and the graph-related (spatial) features from
a range of input data. Those signals include the synthetic time-varying
oscillatory graph signals, resting-state electroencephalogram (EEG) signals,
electrical consumption data of clients within a power network and plant-wide
oscillations from industrial networks. In Section VI, we present results that
compare the proposed method with the only available temporal graph signal
decomposition algorithm (TGSD) method [22] and highlight their crucial
differences.

Figure 1: Time-series plots of input synthetic signal and decomposed
components obtained by applying the TVGMD method.

(a)

(b)

Figure 2: Multi-scale connectivity matrices obtained by applying the TVGMD
method on a) clean input signal (32); b) noisy input signal (32) with the
value of $\eta_{i}$ that correspond to the $SNR=6\mbox{ }dB$ for all $x_{i}$.

###  V-A Analysis of Time-varying Synthetic Graph Signal

The first experiment involves a time-varying synthetic graph signal consisting
of 8 nodes with each containing a combination of oscillations. The goal of
this experiment is to demonstrate that the graph modes obtained from the TVGMD
algorithm exhibit the following properties: i) the graph modes comprise band-
limited time-varying signal oscillations; ii) for each graph mode, the
corresponding adjacency matrix contains information about the connectivity
structure of the mode. Further, the noise robustness of the TVGMD algorithm
will also be assessed on the input data. Owing to the synthetic nature of the
data, the ground truth for each graph mode is readily available that will
facilitate the assessment of the TVGMD method. In this experiment, the
parameter values of the TVGMD method were: $\alpha=200$, $\beta=0.1$,
$\gamma=1$, $\tau=0$ and $K=4$.

The input time-varying graph signal, denoted by
$\mathbf{x}=[x_{1},x_{2},\ldots,x_{8}]$, is defined as follows.

$\displaystyle\centering x_{1}\@add@centering$ $\displaystyle=\cos(2\pi\times
2t)+\cos(2\pi\times 128t)+\eta_{1}$ $\displaystyle x_{2}$
$\displaystyle=-\cos(2\pi\times 24t)+\cos(2\pi\times 48t)+\eta_{2}$
$\displaystyle x_{3}$ $\displaystyle=\cos(2\pi\times 2t)+\cos(2\pi\times
48t)+\eta_{3}$ $\displaystyle x_{4}$ $\displaystyle=\cos(2\pi\times
24t)+\cos(2\pi\times 128t)+\eta_{4}$ $\displaystyle x_{5}$
$\displaystyle=\cos(2\pi\times 2t)+\cos(2\pi\times 24t)+\cos(2\pi\times 48t)$
$\displaystyle\quad+\cos(2\pi\times 128t)+\eta_{5}$ $\displaystyle x_{6}$
$\displaystyle=\cos(2\pi\times 48t)+\cos(2\pi\times 128t)+\eta_{6}$
$\displaystyle x_{7}$ $\displaystyle=\cos(2\pi\times 2t)+\cos(2\pi\times
24t)+\eta_{7}$ $\displaystyle x_{8}$ $\displaystyle=\cos(2\pi\times
2t)+\cos(2\pi\times 24t)-\cos(2\pi\times 48t)+\eta_{8}$ (32)

Note that each of the graph nodes consists of a combination of up to four
oscillatory signals: 2 Hz, 24 Hz, 48 Hz and 128 Hz sinusoid. $\eta_{i}$
denotes the additive Gaussian noise added to the signal associated with the
$i$-th graph node. The oscillations in the graph signal (32) also defines its
connectivity structure at each scale. For instance, as the 2 Hz oscillation is
commonly present within the signals at graph nodes
$[x_{1},x_{3},x_{5},x_{7},x_{8}]$, we expect those nodes to be connected at
the corresponding (low-frequency) scale. At the same scale, the remaining
nodes $[x_{2},x_{4},x_{6}]$ would also be connected as their corresponding
time-series are all close to zero.

The time series plots of the original noiseless graph signal (for
$\eta_{i}=0$) as well as the decomposed graph modes obtained from the TVGMD
method are shown in Fig. 1. Here, the modes $c_{1}$, $c_{2}$, $c_{3}$ and
$c_{4}$ are associated with the oscillations of frequencies 2 Hz, 48 Hz, 24 Hz
and 128 Hz respectively. In each case, the time series plots for all graph
nodes ($N1$-$N8$) are depicted in the Fig. 1 in different colours.

Fig. 2(a) shows the corresponding connectivity structure of each graph mode,
inferred from the corresponding adjacency matrix. Each node is depicted by a
red or a blue marker in the connectivity graphs, with the red (blue) node
indicating that the corresponding graph mode has high (low) signal energy.

Note that there is a complete alignment of information across different scales
in the extracted graph modes as shown in Fig. 1; for instance, only 2 Hz
oscillation is present in $c_{1}$; 48 Hz in $c_{2}$; 24 Hz in $c_{3}$ and 128
Hz in $c_{4}$. This alignment of information is a crucial requirement in many
signal processing applications involving multivariate signals e.g., data
fusion [27] and denoising [28].

Figure 3:  Time series plots of the 12-channel input EEG signal and the
decomposed components obtained by using the TV-GMD algorithm. The x-axis
denotes the number of samples; the sampling frequency was 512 Hz. EEG
amplitudes were scaled in the range of $[-1,1]$ for better illustration.

The analysis of the adjacency matrices of the graph modes in Fig. 2(a) also
reveals the underlying connectivity patterns in the data. For example, the
connectivity pattern of the first graph mode $c_{1}$ suggests that the high-
energy nodes $N1$, $N3$, $N5$, $N7$ and $N8$ are connected, whereas $N2$, $N4$
and $N6$ are close to zero (hence the blue colored nodes) and are thus
connected to each other. This observation is consistent with how the graph
signal has been defined in (32). The same could be concluded for the other
graph modes: in $c_{3}$ for instance, the nodes $N4$, $N5$, $N7$ and $N8$
carry the 24 Hz oscillation signals and are correctly shown as connected,
whereas $N1$, $N3$ and $N6$ are all close to zero and therefore also
connected. An interesting case here is that of node $N2$ that is not connected
to the remaining nodes. The reason is that the term associated with the 24 Hz
sinusoid in $N2$ (or $x_{2}$) (32) has negative sign meaning that it is out-
of-phase with the 24 Hz oscillations at the other nodes.

Finally, the connectivity graphs obtained for the noisy graph signal (32) are
shown in Fig. 2(b). The noise variance at each graph node was set to a level
that resulted in the $SNR=6\mbox{ }dB$ for each of the input signal $x_{i}$.
Comparing these connectivity graphs with those obtained for the clean signal
in Fig. 2(a), it is clear that the addition of noise affects the output of the
TVGMD method and hence alters the connectivity patterns. That said, the
primary data patterns at multiple scales remain intact. Specifically, for
$c_{2}$, we observe a partially-connected graph between $N2$, $N3$, $N5$ and
$N6$ instead of a fully connected one. Similar observations can be made
regarding the connectivity patterns of $c_{4}$. While this result gives some
insight into the performance of the TVGMD in the presence of noise, a detailed
study on the topic could be an avenue for future work.

Figure 4:  Connectivity graphs of the EEG data at multiple scales obtained
from TVGMD. No connectivity exists among the nodes at $c_{1}$, whereas strong
connectivity exists among the nodes representing the occipital and the
parietal regions within $c_{2}$. This makes sense as $c_{2}$ carries alpha
rhythms which are typically dominant in the occipital and parietal brain
regions. Figure 5:  Multi-mode connectivity graphs of the electricity
consumption time-varying graph signal, obtained via TVGMD.

###  V-B Analysis of Alpha Waves in Electroencephalogram (EEG)

This experiment will demonstrate the potential of the TVGMD method in the
analysis of biomedical electroencephalogram (EEG) signals. EEG is a recording
of brain’s electrical activity, obtained from multiple sensors (electrodes)
attached to the scalp. Typically, each electrode records a time-series of
electrical signal produced by brain. This way, multichannel EEG signals can be
represented by time-varying graph signals with each electrode representing a
graph node.

EEG signals are known to exhibit oscillatory behaviour in specific frequency
bands e.g., theta (4-8Hz), alpha (8-12 Hz), beta (13-30 Hz) and gamma (30-150
Hz) bands, each linked with different underlying neural processes. Therefore,
a crucial task in numerous applications involving EEG signals is to decompose
neural oscillations into specific bands. Further, it is of interest to reveal
the functional connectivity structures in brain at multiple frequency bands.
Here, we demonstrate the potential of the TVGMD method to concurrently provide
both pieces of information i.e., time-varying modes and their connectivity
structures.

Here, we apply the TVGMD method to the publicly available EEG data that was
obtained in a resting-state eyes-open/closed experimental protocol. The data
was recorded in the GIPSA-lab, Grenoble, France, in 2017. A research-grade
amplifier (g.USBamp, g.tec, Schiedlberg, Austria) along with an EC20 cap
comprising 16 electrodes, placed according to the 10-20 international system,
was used to collect the data. Further details about the data acquisition
process and the related experiment can be found in [29]. The TVGMD parameter
values used in this experiment were: $\alpha=1000$, $\beta=0.5$, $\gamma=1$,
$\tau=0.1$ and $K=6$.

The results of applying the TVGMD algorithm on the 12-channel EEG signal is
shown. Particularly, the input and the time-series associated with the
obtained graph modes are shown in the Fig. 3 while the corresponding
functional connectivity structures are shown in the Fig. 4. Note that the time
plots of the first mode ($c_{1}$) mainly contains the low-frequency trend of
the signal. A look at the corresponding functional connectivity map of $c_{1}$
in the lower figure reveals no prominent connections among different brain
regions as expected. The power spectrum of the graph modes corresponding to
each node are shown in the inset plot with each connectivity graph.

The time series of the second graph mode, $c_{2}$, comprises oscillations
related to the alpha band which are of particular interest here since the EEG
data was obtained in the eyes-closed state, where alpha rhythms are typically
prominent in the occipital and the parietal regions. This can be verified in
the corresponding connectivity map of $c_{2}$ that not only shows consistently
higher signal power in the occipital and the parietal regions but also a
strong functional connectivity within that region. Note that in the
corresponding Fourier spectra of the graph mode (shown in the inset of the
connectivity graph), there is a peak at the 8-12 Hz frequency range which
implies that we are indeed operating in the alpha frequency band.

For the remaining couple of graph modes, $c_{3}$ and $c_{4}$, the time plots
and the connectivity patterns are also shown. For these modes, we observe
mostly low-energy signals as well as random connectivity patterns in the
different brain regions.

Figure 6: Time-series plots of the multiple graph modes of the electricity
consumption data at five selected node.

###  V-C Electric Consumption Data Analysis in Power Networks

The next time-varying graph data set used in our experiment is the electricity
consumption in kWh, recorded every 15 minutes from 2012 to 2014 in Portugal
for 206 clients [30]. The data was transformed to show hourly consumption of
electricity and 480 hours (20 days) of data was analyzed using the TVGMD
method; the results are shown in Fig. 5 and Fig. 6. The following TVGMD
parameter values were used in this experiment: $\alpha=2000$, $\beta=0.25$,
$\gamma=1$, $\tau=0$ and $K=10$.

Fig. 5 shows the connectivity structures of the first 6 graph modes obtained
from the TVGMD method. Note that 5 selected clients (nodes) have been
specifically highlighted in the connectivity graphs: N7, N19, N78, N105 and
N162. The electricity consumption time-series plots of the graph modes for
these 5 clients (nodes) along with the overall consumption are shown in Fig.
6. These nodes have been carefully selected out of the total 206 nodes to show
a variety of consumption patterns, ranging from a high consumption (at Nodes
N78 and N105), mid-range consumption (at the Node N162) to a very low
consumption (at Nodes N7 and N19).

By analyzing the graph modes in time domain, it was noticed that the
electricity consumption followed oscillatory patterns at different time scales
that were extracted accurately through the obtained modes. Specifically,
$c_{5}-c_{1}$ provided information about the 6-hourly, 8-hourly, 12-hourly,
daily and 1-weekly electricity consumption patterns respectively. The
adjacency graphs gave further information regarding the connectivity of
different nodes (clients) at different time scales in terms of their
electricity usage: for example, the $c_{1}$ and $c_{2}$ time plots show that
the nodes 78 and 105 follow similar consumption patterns on weekly and daily
(24 hourly) basis which can be verified through their strong links in the
$c_{1}$ and $c_{2}$ connectivity graphs. The nodes N7 and N19 also exhibited
strong connections (in terms of electricity usage) across a range of different
time scales that could be verified from the respective connectivity graphs
(though those connections remain mostly hidden in our illustration due to
their closeness and a very high density of the other nodes).

###  V-D Detection of Plant-wide Oscillations in Flotation Circuits

Flotation is used to separate useful mineral particles from gangue particles
in a concentrator process [31]. This is accomplished by imposing appropriate
conditions that enable hydrophobic materials to attach to air bubbles and rise
to the top of the flotation cell. Here, we consider data from a flotation cell
that consists of two parallel banks with seven cells in series in each bank;
the resulting signal can be modelled as a time-varying graph signal with 14
nodes. The process flow diagram of such a flotation cell is shown in Fig. 1 in
[31] along with the relevant description. The time-series data corresponding
to the cell levels and outflows in the flotation circuit are of interest for
analyzing the flotation plant process. Particularly, it is important to
automatically detect oscillations in the plant data as well as diagnosing
their root cause [32].

Figure 7:  Time-series plots of the decomposed graph modes obtained by
applying the TVGMD algorithm on the level signals in a flotation circuit.
Figure 8: Connectivity graphs corresponding to the graph modes from Fig. 7. A
natural grouping of the nodes (flotation cells) is evident from the $c_{2}$
and the $c_{3}$ connectivity graphs.

The graph modes obtained from applying the TVGMD algorithm on input data are
shown in Fig. 7 and Fig. 8. The TVGMD parameter values used in this experiment
were: $\alpha=1000$, $\beta=0.1$, $\gamma=1$, $\tau=0$ and $K=4$. The input
data comprise the level time-series, $L_{1}-L_{14}$, recorded from the 14
cells of the flotation circuit and is shown in the top row of Fig. 7. The time
plots of the decomposed graph modes clearly show the presence of two separate
frequencies in the two flotation banks. Note that the lower frequency
oscillation is prominent in the level time series of the first seven cells (in
$c_{2}$), whereas the higher frequency oscillation is predominantly present in
the remaining cells ($L_{8}-L_{14}$), in $c_{3}$. There is also a clear
alignment of information present in the graph modes e.g., $c_{2}$ ($c_{3}$)
contain only the specific low (high) frequency oscillation in all the cells or
nodes of the graph signal. This mode-alignment is a crucial requirement for
the processing of multivariate signals in many application areas [27].

The corresponding connectivity graphs of the extracted graph modes are shown
in Fig. 8. The connectivity graphs corresponding to $c_{2}$ and $c_{3}$ are of
particular interest here. The $c_{2}$ connectivity graph shows that the level
signals $L_{2}$-$L_{7}$ (from the first bank) are high-energy, as evident from
the node colors, and are correlated with each other. The correlations are weak
though (as evident from the thin connectivity links) owing to the inherent
phase difference between the signals from different nodes. On the other hand,
the connectivity graphs corresponding to $c_{3}$ show that the level signals
$L_{9}$-$L_{14}$ (from the second bank) have higher amplitudes and are
correlated (albeit weakly) to each other. The connectivity graphs can
therefore help infer relationships between measured variables at multiple
scales. This could be used to devise grouping strategies for the diagnosis of
root cause of oscillations in the flotation circuit, facilitating their fault
diagnosis.

##  VI Discussions

The TVGMD method leverages the power of graph-theoretic tools to obtain
inherent oscillatory signal components along with their network connectivity
structures – termed here as the graph modes. A related approach, TGSD [22],
also uses graph-related concepts to extract oscillatory signal components.
Here we compare the results of the two methods on the electricity consumption
data; the data set is explained in Section V-C and is plotted as the top left
subplot of the Fig. 6. The oscillatory modes from the TVGMD and the TGSD
methods are shown in Fig. 6 and Fig. 9 respectively. The TGSD results shown
here are based on its best-tuned parameters so that the obtained modes are
physically meaningful. As noted earlier, each (graph) mode from the TVGMD
method corresponds to the electricity consumption pattern at a specific time
scale e.g., $c_{5}-c_{1}$ are associated with 6-hourly, 8-hourly, 12-hourly,
daily and 1-weekly consumption patterns respectively. In contrast, the TGSD
modes could not be associated to specific time scales and there is a clear
evidence of multiple oscillations appearing in a single mode e.g., in $c_{1}$
for the $N78$ data. One can also observe similar information spanning multiple
modes e.g., across $c_{3}-c_{5}$. For these reasons, it is difficult to
extract meaningful information about users’ electricity consumption patterns
from the TGSD modes. The improved performance of TVGMD over TGSD could be
attributed to the data-driven nature of the TVGMD i.e., no fixed dictionary
atoms (basis functions) are used in the decomposition. More importantly, the
TVGMD method extracts the multi-scale network connectivity patterns of the
data concurrently with the oscillatory modes. TGSD and no other method, to our
knowledge, is able to accomplish that.

Figure 9: TGSD output from the electricity consumption data from Section V-C.

TVGMD method has been generally found to be reasonably robust to the changes
in its parameters, except for the number of graph modes $K$. In Table I,
suggested range of values for each parameter of TVGMD has been provided. The
parameter $\alpha$ dictates the bandwidth of the graph modes: the lower its
value, the higher the bandwidth of the modes. In our experiments, the range of
values $\alpha=500-2000$ provided good results, depending on the whether the
low-bandwidth or high-bandwidth graph modes are needed. Lower values of the
parameter $\beta$ i.e., in the range of $0.1-1$, and $\gamma$, within the
range of $0.5-1.0$, were found to produce optimal graph modes in our
experiments. The nonzero value of the parameter $\tau$ enables the term
involving Lagrangian multiplier in the proposed model, thus enforcing exact
reconstruction of input signal via graph modes. The lower values of $\tau$
(including $\tau=0$), therefore, are useful in the case of noisy input signal
where its exact reconstruction is not desirable. Finally, the choice of the
number of graph modes $K$ was mainly dictated by the input data at hand and
the performance of the algorithm depended heavily on the correct choice of
$K$. For instance, in the case of synthetic graph signal (32) for which $K=4$
is known by design, choosing any other value of $K$ resulted in significant
graph mode-mixing (multiple oscillatory components in a single mode). This
issue represents one of the limitations of the TVGMD method as the knowledge
of correct value of $K$ may not be known a priori.

While the convergence of the TVGMD method has not been established and is not
in the scope of this work, the TVGMD method converged in all the experiments
shown in this work as well as in other simulations involving a wide range of
signals.

TABLE I: User-defined parameter values within TVGMD Parameters | $\alpha$ | $\beta$ | $\gamma$ | $\tau$ | $K$
---|---|---|---|---|---
TVGMD | 500-2000 | 0.1-1 | 0.5-1 | 0-1 | input dependent

##  VII Conclusion

In this paper, a new method for extracting the inherent modes of time-varying
graph signals has been proposed. The extracted graph modes are defined by the
following useful pieces of information: i) the graph modes comprise inherent
oscillatory signals present in the data at multiple nodes; ii) there is an
alignment of information provided by the graph modes at each node (vertex);
iii) the graph modes convey information about the underlying connectivity
networks of the time-varying graph signals at multiple frequency scales. The
method is the first-of-its-kind to provide the above information about the
temporal graph signals concurrently by using a fully data-driven approach. By
data-driven approach, it is meant that no predefined basis functions (either
in the time- or graph-domains) have been used to obtain the graph modes.

At the heart of the proposed method is a robust variational optimization
formulation that includes multiple requirements (or constraints) for the
extraction of finite number of graph modes from an input graph signal. Those
requirements relate to the properties of the graph modes both in the time-
domain (e.g., limited bandwidth, complete signal reconstruction) and along the
graph geodesics (e.g., smoothness of the graph modes along their respective
adjacency matrices). We have proposed a method to solve the resulting
optimization problem by using the alternative direction method of multipliers
(ADMM) and the primal-dual technique. The proposed method outputs graph modes
that comprise the time-series associated with each node at multiple frequency
scales and their corresponding network connectivity structure.

To demonstrate the power and the potential of the proposed method in different
domains, results from multiple experiments have been included in the paper.
Particularly, the proposed algorithm has been shown to isolate alpha rhythms
from the resting-state EEG signals and to correctly localize their sources in
the occipital and parietal regions in the brain (See Fig. 4). Further, the
proposed method was used to obtain the connectivity patterns of users in terms
of their electricity consumption at different time scales (Fig. 5). Finally,
the connectivity graphs of the flotation plant data set (Fig. 8) could be used
to define a grouping strategy of the plant cells, facilitating the fault
diagnosis in flotation plants.

Future work could consider extending this method to obtain time-varying
adjacency (connectivity) matrices. Currently, the method obtains static
connectivity matrices based on the graph signal evolution for the entire
signal duration. This is clearly a limitation of the proposed method since
graph signals can exhibit non-stationarity in terms of time-varying
connectivity patterns. Other avenues for future work may include: i) the
development of computationally efficient approaches that target the extraction
of a single graph mode at a time; ii) new adaptive approaches that
automatically set the optimal number of graph modes $K$.

##  VIII Appendix

To solve the convex optimization problem given in (28), a primal-dual
splitting approach is employed that jointly solves the primal and the
associated dual optimization problems instead of focusing exclusively on
either one. These algorithms are discussed in detail in [26]. Among these
class of algorithms, the one used in this work is the forward-backward-forward
(FBF) based primal-dual algorithm that combines the gradient descent step(s)
with the computation step(s) involving the proximity operator. The steps of
that algorithm are listed in Algorithm 6 in [26]. In what follows, we list the
steps to solve (28) by using the FBF based primal-dual method.

Algorithm 3 Primal-Dual Algorithm to Solve (28)

Input:
$z,\beta,\gamma,w_{1}^{(k)},d_{1}\in\mathbb{R}_{+}^{m},\delta,\epsilon,i\leftarrow
0$

repeat

$i\leftarrow i+1$

$y_{i+1}\leftarrow w_{i}-\delta\times(2\gamma w_{i}+Q^{T}d_{i})$

$\mkern 1.5mu\overline{\mkern-1.5muy\mkern-1.5mu}\mkern 1.5mu_{i+1}\leftarrow
d_{i}+\delta\times(Qw_{i})$

$p_{i+1}\leftarrow\max(0,y_{i+1}-2\beta\delta z)$

$\mkern 1.5mu\overline{\mkern-1.5mup\mkern-1.5mu}\mkern
1.5mu_{i+1}\leftarrow\frac{(\mkern
1.5mu\overline{\mkern-1.5muy\mkern-1.5mu}\mkern 1.5mu_{i+1}-\sqrt{\mkern
1.5mu\overline{\mkern-1.5muy\mkern-1.5mu}\mkern 1.5mu_{i+1}^{2}+4\delta})}{2}$

$q_{i+1}\leftarrow p_{i+1}-\delta\times(2\gamma p_{i+1}+Q^{T}\mkern
1.5mu\overline{\mkern-1.5mup\mkern-1.5mu}\mkern 1.5mu_{i+1})$

$\mkern 1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern
1.5mu_{i+1}\leftarrow\mkern 1.5mu\overline{\mkern-1.5mup\mkern-1.5mu}\mkern
1.5mu_{i+1}+\delta\times(Qp_{i+1})$

$(w_{i+1},d_{i+1})=(w_{i}-y_{i+1}+q_{i+1},w_{i}-\mkern
1.5mu\overline{\mkern-1.5muy\mkern-1.5mu}\mkern 1.5mu_{i+1}+\mkern
1.5mu\overline{\mkern-1.5muq\mkern-1.5mu}\mkern 1.5mu_{i+1})$

until Convergence:
$\frac{\|w_{i+1}-w_{i}\|_{2}}{\|w_{i}\|_{2}}<\epsilon\And\frac{\|d_{i+1}-d_{i}\|_{2}}{\|d_{i}\|_{2}}<\epsilon$

## References

  * [1] A. Ortega, P. Frossard, J. Kovačević, J. M. F. Moura, and P. Vandergheynst, “Graph signal processing: Overview, challenges, and applications,” _Proc. IEEE_ , vol. 106, no. 5, pp. 808–828, 2018.
  * [2] W. Huang, T. A. W. Bolton, J. D. Medaglia, D. S. Bassett, A. Ribeiro, and D. Van De Ville, “A graph signal processing perspective on functional brain imaging,” _Proc. IEEE_ , vol. 106, no. 5, pp. 868–885, 2018.
  * [3] B. Zhao, K. He, L. Stankovic, and V. Stankovic, “Improving event-based non-intrusive load monitoring using graph signal processing,” _IEEE Access_ , vol. 6, pp. 53 944–53 959, 2018.
  * [4] I. Jabłoński, “Graph signal processing in applications to sensor networks, smart grids, and smart cities,” _IEEE Sensors J._ , vol. 17, no. 23, pp. 7659–7666, 2017.
  * [5] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on graphs,” _IEEE Trans. Signal Process._ , vol. 61, no. 7, pp. 1644–1656, 2013.
  * [6] N. Tremblay and P. Borgnat, “Subgraph-based filterbanks for graph signals,” _IEEE Trans. Signal Process._ , vol. 64, no. 15, pp. 3827–3840, 2016.
  * [7] W. Xu, E. Mallada, and A. Tang, “Compressive sensing over graphs,” in _2011 Proc. IEEE INFOCOM_ , 2011, pp. 2087–2095.
  * [8] N. Perraudin and P. Vandergheynst, “Stationary signal processing on graphs,” _IEEE Trans. Signal Process._ , vol. 65, no. 13, pp. 3462–3477, 2017.
  * [9] L. Stanković, J. Lerga, D. Mandic, M. Brajović, C. Richard, and M. Daković, “From time–frequency to vertex–frequency and back,” _Mathematics_ , vol. 9, no. 12, 2021.
  * [10] 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,” _IEEE Trans. Signal Process._ , vol. 66, no. 3, pp. 817–829, 2018.
  * [11] J. H. Giraldo, A. Mahmood, B. Garcia-Garcia, D. Thanou, and T. Bouwmans, “Reconstruction of time-varying graph signals via sobolev smoothness,” _IEEE Trans. Signal Inf. Process. Netw._ , vol. 8, pp. 201–214, 2022.
  * [12] J. Jiang, H. Feng, D. B. Tay, and S. Xu, “Theory and design of joint time-vertex nonsubsampled filter banks,” _IEEE Trans. Signal Process._ , vol. 69, pp. 1968–1982, 2021.
  * [13] D. B. Tay and J. Jiang, “Time-varying graph signal denoising via median filters,” _IEEE Trans. Circuits Syst. II_ , vol. 68, no. 3, pp. 1053–1057, 2021.
  * [14] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, S. H. H., Q. Zheng, N.-C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis,” _Proc. R. Soc. A: Math. Phys. Eng. Sci._ , vol. 454, no. 1971, pp. 903–995, 1998.
  * [15] K. Dragomiretskiy and D. Zosso, “Variational mode decomposition,” _IEEE Trans. Signal Process._ , vol. 62, no. 3, pp. 531–544, 2014.
  * [16] I. Daubechies, J. Lu, and H.-T. Wu, “Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool,” _Appl. Comput. Harmon. Anal._ , vol. 30, no. 2, pp. 243–261, 2011.
  * [17] J. Harmouche, D. Fourer, F. Auger, P. Borgnat, and P. Flandrin, “The sliding singular spectrum analysis: A data-driven nonstationary signal decomposition tool,” _IEEE Trans. Signal Process._ , vol. 66, no. 1, pp. 251–263, 2018\.
  * [18] N. Rehman and D. Mandic, “Multivariate empirical mode decomposition,” _Proc. R. Soc. A: Math. Phys. Eng. Sci._ , vol. 466, no. 2117, pp. 1291–1302, 2010.
  * [19] N. Rehman and H. Aftab, “Multivariate variational mode decomposition,” _IEEE Trans. Signal Process._ , vol. 67, no. 23, pp. 6039–6052, 2019.
  * [20] A. Cicone and E. Pellegrino, “Multivariate fast iterative filtering for the decomposition of nonstationary signals,” _IEEE Trans. Signal Process._ , vol. 70, pp. 1521–1531, 2022.
  * [21] T. Eriksen and N. Rehman, “Data-driven signal decomposition approaches: A comparative analysis,” 2022. [Online]. Available: https://arxiv.org/abs/2208.10874
  * [22] M. J. McNeil, L. Zhang, and P. Bogdanov, “Temporal graph signal decomposition,” in _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_. ACM, aug 2021.
  * [23] Y. Yankelevsky and M. Elad, “Finding gems: Multi-scale dictionaries for high-dimensional graph signals,” _IEEE Trans. Signal Process._ , vol. 67, no. 7, pp. 1889–1901, 2019.
  * [24] V. Kalofolias, “How to learn a graph from smooth signals,” in _Proc. IASTATS Int. Conf. Artif Intell. Stat._ , ser. Proceedings of Machine Learning Research, A. Gretton and C. C. Robert, Eds., vol. 51. Cadiz, Spain: PMLR, 09–11 May 2016, pp. 920–929.
  * [25] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Found. Trends Mach. Learn._ , vol. 3, no. 1, pp. 1–122, 2011\.
  * [26] N. Komodakis and J.-C. Pesquet, “Playing with duality: An overview of recent primal-dual approaches for solving large-scale optimization problems,” _IEEE Signal Process. Mag._ , vol. 32, no. 6, pp. 31–54, 2015.
  * [27] N. Rehman, S. Ehsan, S. M. U. Abdullah, M. J. Akhtar, D. P. Mandic, and K. D. McDonald-Maier, “Multi-scale pixel-based image fusion using multivariate empirical mode decomposition,” _Sensors_ , vol. 15, no. 5, pp. 10 923–10 947, 2015.
  * [28] N. Rehman, B. Khan, and K. Naveed, “Data-driven multivariate signal denoising using mahalanobis distance,” _IEEE Signal Process. Lett._ , vol. 26, no. 9, pp. 1408–1412, 2019.
  * [29] G. Cattan, P. L. C. Rodrigues, and M. Congedo, “Eeg alpha waves dataset,” Dec. 2018. [Online]. Available: https://doi.org/10.5281/zenodo.2605110
  * [30] D. Dua and C. Graff, “UCI machine learning repository,” 2017. [Online]. Available: http://archive.ics.uci.edu/ml
  * [31] B. Lindner, M. Chioua, J. Groenewald, L. Auret, and M. Bauer, “Diagnosis of oscillations in an industrial mineral process using transfer entropy and nonlinearity index,” _IFAC-PapersOnLine_ , vol. 51, no. 24, pp. 1409–1416, 2018.
  * [32] Q. Chen, X. Lang, S. Lu, N. Rehman, L. Xie, and H. Su, “Detection and root cause analysis of multiple plant-wide oscillations using multivariate nonlinear chirp mode decomposition and multivariate granger causality,” _Comput. Chem. Eng._ , vol. 147, p. 107231, 2021.