# Lossless digraph signal processing via polar decomposition

Feng Ji The author is with the School of Electrical and Electronic
Engineering, Nanyang Technological University, 639798, Singapore (e-mail:
jifeng@ntu.edu.sg).

###### Abstract

In this paper, we present a signal processing framework for directed graphs.
Unlike undirected graphs, a graph shift operator such as the adjacency matrix
associated with a directed graph usually does not admit an orthogonal
eigenbasis. This makes it challenging to define the Fourier transform. Our
methodology leverages the polar decomposition to define two distinct
eigendecompositions, each associated with different matrices derived from this
decomposition. We propose to extend the frequency domain and introduce a
Fourier transform that jointly encodes the spectral response of a signal for
the two eigenbases from the polar decomposition. This allows us to define
convolution following a standard routine. Our approach has two features: it is
lossless as the shift operator can be fully recovered from factors of the
polar decomposition. Moreover, it subsumes the traditional graph signal
processing if the graph is directed. We present numerical results to show how
the framework can be applied.

###### Index Terms:

Directed graph, graph signal processing, polar decomposition

##  I Introduction

The field of graph signal processing (GSP) is experiencing significant growth
and focuses on the analysis of signals defined on graphs [1, 2, 3, 4, 5, 6, 7,
8]. Numerous real-world phenomena, including social networks, transportation
systems, and sensor networks, can be naturally represented as graphs. In GSP,
the key idea revolves around the vector space of graph signals, where each
node in a given graph is assigned a numerical value by a graph signal. As
signals on a graph form a vector space, GSP employs linear transformations,
such as the graph Fourier transform (GFT) and graph filters, to examine and
explore relationships among graph signals.

On the other hand, despite a growing need for applications that involve
directed graphs (digraphs) [9, 10], the intrinsic asymmetry of their edge
connections presents challenges in defining GFTs. First, the basic graph
operators become non-diagonalizable for a large number of digraphs. This makes
it challenging to define key concepts such as the frequency domain and
convolution. Second, “ordering" the graph Fourier basis, so that we can
interpret them as low or high frequencies is difficult. This is because the
eigenvalues, which can be used for the ordering, can be complex numbers.

As a consequence, a variety of approaches have been proposed to address these
challenges. For example, studies directly using the adjacency matrix as the
graph operator have suggested utilizing the Jordan decomposition [11, 12] or
the singular value decomposition [13] of the adjacency matrix for constructing
GFT. Symmetrization of the adjacency matrix has also been studied [14]. [15]
proposes solving an optimization based on graph directed variation to obtain
an orthonormal basis for Fourier transform. There are studies that define the
graph operator as decomposed matrices of the adjacency matrix, producing
diagonalizable matrices [16, 17]. In [18], the authors propose to directly
modify the graph Laplacian so that the resulting one is Hermitian and hence
diagonalizable. The recent work [19] proposes to modify the graph structure so
that the corresponding graph operator on the resulting graph becomes
diagonalizable. However, these approaches have their respective shortcomings.
For example, the Jordan decomposition is not orthogonal. Hence, signal energy
is not preserved after the Fourier transformation, and the basis vectors have
non-trivial interactions among themselves. The approaches of optimization and
topology perturbation are usually not lossless. The original directed graph
topology cannot be faithfully recovered from the basis of the Fourier
transform.

In this paper, we propose a lossless digraph signal process scheme. There are
two main ingredients of our approach. For a graph shift operator (GSO)
$\mathbf{S}$, which may not be symmetric, we consider the polar decomposition
of $\mathbf{S}$. This means that $\mathbf{S}$ is the product of a unitary
matrix $\mathbf{U}$ and symmetric matrix $\mathbf{P}$, both are orthogonally
diagonalizable. Then, we introduce the Fourier transform that utilizes the
orthonormal basis of $\mathbf{U}$ and $\mathbf{P}$ jointly. Unlike the
traditional GSP, the frequency domain is $n^{2}$ dimensional instead of $n$
dimensional, which allows us to package more information and thus the
framework is lossless. Moreover, we also demonstrate that our framework
subsumes the traditional GSP if $\mathbf{S}$ is symmetric, and hence we have a
strict generalization of GSP.

The rest of the paper is organized as follows. In Section II, we review the
traditional GSP and matrix polar decomposition, which is the fundamental tool
of this paper. In Section III, we introduce the Fourier transform. The main
idea is to extend the frequency domain from $\mathbb{R}^{n}$ to
$\mathbb{C}^{n^{2}}$. Once we have properly introduced the frequency domain,
we define the notion of convolution in Section IV following a standard
routine. In Section V, we study how matrix perturbation affects the framework.
We discuss miscellaneous topics in Section VI such as how to resolve the issue
of non-uniqueness of polar decomposition and eigendecomposition. We present
numerical results in Section VII and conclude in Section VIII.

##  II Preliminaries

In this section, we review the basics of graph signal processing (GSP) [1] and
polar decomposition [20]. In the paper, we use the term “GSP” exclusively for
the traditional graph signal processing on undirected graphs.

Let $G=(V,E)$ be a graph of size $n$. If $G$ is undirected, GSP requires a
symmetric graph shift operator (GSO) $\mathbf{S}$. Common choices of
$\mathbf{S}$ include its adjacency matrix $\mathbf{A}$ and the Laplacian
matrix $\mathbf{L}=\mathbf{A}-\mathbf{D}$, where $\mathbf{D}$ is the degree
matrix. By the spectral theorem, the symmetric matrix $\mathbf{S}$ admits an
eigendecomposition
$\mathbf{S}=\mathbf{V}\bm{\Lambda}\mathbf{V}^{{}^{\intercal}}$, where
$\bm{\Lambda}$ is the diagonal matrix consisting of eigenvalues of
$\mathbf{S}$ and the columns of $\mathbf{V}$ are the corresponding
eigenvectors. Thus $\mathbf{V}$ is an orthogonal matrix. Suppose
$V=\\{v_{1},\ldots,v_{n}\\}$ has a fixed ordering. A graph signal
$\mathbf{f}=(x_{i})_{1\leq i\leq n}$ is a vector in $\mathbb{R}^{n}$, where
$x_{i}$ is the signal assigned to $v_{i}$. The graph Fourier transform (GFT)
of $\mathbf{x}$ w.r.t. $\mathbf{S}$ is

$\displaystyle\widehat{\mathbf{f}}=\mathbf{V}^{{}^{\intercal}}\mathbf{f}\in\mathbb{R}^{n}.$
(1)

The frequency domain is indexed by the eigenvalues of $\mathbf{S}$.

However, if $G$ is directed, then we have the challenge that
$\mathbf{S}=\mathbf{A}$ or $\mathbf{L}$ may not be symmetric. We usually do
not have an eigendecomposition and thus important signal processing notions
such as GFT (1) cannot be properly defined. We propose to resolve this issue
using matrix polar decomposition [20, 21], which is reviewed next.

Recall that for any $n\times n$ matrix $\mathbf{S}$, its polar decomposition
is $\mathbf{S}=\mathbf{U}\mathbf{P}$, where $\mathbf{U}$ is an orthogonal
matrix and $\mathbf{P}$ is symmetric and positive semi-definite. Moreover, if
$\mathbf{S}$ is invertible, then the decomposition is unique, and $\mathbf{P}$
is positive definite. Explicitly, let
$\mathbf{S}=\mathbf{V}\bm{\Sigma}\mathbf{W}^{{}^{\intercal}}$ be the singular
value decomposition of $\mathbf{S}$ [22], then
$\mathbf{U}=\mathbf{V}\mathbf{W}^{{}^{\intercal}}$ and
$\mathbf{P}=\mathbf{W}\bm{\Sigma}\mathbf{W}^{{}^{\intercal}}$.

Both $\mathbf{U}$ and $\mathbf{P}$ are normal and thus they admit (complex)
eigendecompositions, i.e.,
$\mathbf{U}=\mathbf{V}_{1}\bm{\Lambda}_{1}\mathbf{V}_{1}^{H}$ and
$\mathbf{P}=\mathbf{V}_{2}\bm{\Lambda}_{2}\mathbf{V}_{2}^{H}$, where
$\mathbf{V}_{1},\mathbf{V}_{2}$ are unitary,
$\bm{\Lambda}_{1},\bm{\Lambda}_{2}$ are diagonal and H is the adjoint operator
(Hermitian transpose). Geometrically, $\mathbf{P}$ performs scaling, while
$\mathbf{U}$ performs a rotation/reflection. By convention, we increasingly
order the eigenvalues of $\bm{\Lambda}_{2}$ according to the absolute value
and those of $\bm{\Lambda}_{1}$ according to the complex phase. We shall use
$\mathbf{V}_{1},\mathbf{V}_{2}$ to define the Fourier transform in the next
section.

##  III The Fourier transform

In this section, we define the Fourier transform. We use the same notations as
in Section II for the polar decomposition. To motivate, consider a (complex)
graph signal $\mathbf{f}$. We want to study its shifted signal
$\mathbf{g}=\mathbf{S}(\mathbf{f})=\mathbf{U}\mathbf{P}(\mathbf{f})$. Its
magnitude $|\mathbf{g}|$ is the same as $|\mathbf{P}(\mathbf{f})|$. Therefore,
taking ordinary GFT (1) w.r.t. $\mathbf{P}$ gives us “smoothness” information
of $\mathbf{f}$. However, rotational information is missing if we ignore
$\mathbf{U}$. To incorporate contributions from both $\mathbf{P}$ and
$\mathbf{U}$, we follow the idea of [8] and propose to expand the spectral
domain from $\mathbb{R}^{n}$ to
$\mathbb{R}^{n}\otimes\mathbb{C}^{n}\subset\mathbb{C}^{n}\otimes\mathbb{C}^{n}\cong\mathbb{C}^{n^{2}}$.
Following the ordering of the matrix multiplication, the proposed Fourier
transform should encode how $\mathbf{f}$ respond w.r.t. $\mathbf{P}$ and how
eigenvectors of $\mathbf{P}$ should respond w.r.t. $\mathbf{U}$ (cf. Section
VI-C below).

Specifically, we denote the columns of $\mathbf{V}_{1}$ by
$\\{\mathbf{v}_{1,1},\ldots,\mathbf{v}_{1,n}\\}$ and of $\mathbf{V}_{2}$ by
$\\{\mathbf{v}_{2,1},\ldots,\mathbf{v}_{2,n}\\}$. They are the eigenvectors of
$\mathbf{V}_{1}$ and $\mathbf{V}_{2}$ respectively. For the signal
$\mathbf{f}$, we may express it as

$\displaystyle\mathbf{f}=\sum_{1\leq i\leq
n}\widehat{\mathbf{f}}_{i}\mathbf{v}_{2,i},$ (2)

where
$\widehat{\mathbf{f}}_{i}=\langle\mathbf{v}_{2,i},\mathbf{f}\rangle=\mathbf{v}_{2,i}^{H}\mathbf{f}$
is the GFT of $\mathbf{f}$ w.r.t. $\mathbf{P}$ and $\langle\cdot,\cdot\rangle$
is the complex inner product.

On the other hand, for each $\mathbf{v}_{2,i}$, we may inspect its “angle”
w.r.t. to each direction of rotation $\mathbf{v}_{1,j}$ as
$p_{i,j}=\langle\mathbf{v}_{1,j},\mathbf{v}_{2,i}\rangle$. Therefore,

$\displaystyle\mathbf{v}_{2,i}=\sum_{1\leq j\leq n}p_{i,j}\mathbf{v}_{1,j}.$
(3)

Combining (2) and (3), we have

$\displaystyle\mathbf{f}=\sum_{1\leq j\leq n}\Big{(}\sum_{1\leq i\leq
n}\widehat{\mathbf{f}}_{i}p_{i,j}\Big{)}\mathbf{v}_{1,j}.$

Therefore, we define the Fourier transform by considering contributions from
$\mathbf{V}_{1}$ and $\mathbf{V}_{2}$ jointly.

######  Definition 1.

The _Fourier transform_ w.r.t. $\mathbf{S}$ is
$\mathcal{F}_{\mathbf{S}}:\mathbb{R}^{n}\to\mathbb{C}^{n^{2}}$ is defined as
follows: for a signal $\mathbf{f}$, then

$\displaystyle\mathcal{F}_{\mathbf{S}}(\mathbf{f})_{i,j}=\widehat{\mathbf{f}}_{i}p_{i,j}=\langle\mathbf{v}_{2,i},\mathbf{f}\rangle\langle\mathbf{v}_{1,j},\mathbf{v}_{2,i}\rangle.$
(4)

For simplicity, we usually denote $\mathcal{F}_{\mathbf{S}}(\mathbf{f})$ by
$\mathcal{F}(\mathbf{f})$ or $\widetilde{\mathbf{f}}$ if $\mathbf{S}$ is clear
from the context.

The _inverse transform_
$\mathcal{I}_{\mathbf{S}}:\mathbb{C}^{n^{2}}\to\mathbb{C}^{n}$ (or simply
$\mathcal{I}$) is

$\displaystyle(a_{i,j})_{1\leq i,j\leq n}\mapsto\sum_{1\leq j\leq
n}\sum_{1\leq i\leq n}a_{i,j}\mathbf{v}_{1,j}.$

To give an explicit matrix formula for the Fourier transform, let $\mathbf{g}$
be any (column) vector, and use $D(\mathbf{g})$ to denote the diagonal matrix
with diagonal $\mathbf{g}$. Moreover, use $\mathbf{c}$ to denote the constant
(column) vector with each entry $1$. Then we have the following observations.

######  Lemma 1.

  1. (a)

$\mathcal{F}(\mathbf{f})=\mathbf{V}_{1}^{H}\mathbf{V}_{2}D(\mathbf{V}_{2}^{H}\mathbf{f})$
and $\mathcal{I}(\mathbf{M})=\mathbf{V}_{1}\mathbf{M}\mathbf{c}$ for
$\mathbf{f}\in\mathbb{R}^{n},\mathbf{M}\in\mathbb{C}^{n^{2}}$.

  2. (b)

Parseval’s identity:
${\left\lVert{\mathbf{f}}\right\rVert}={\left\lVert{\mathcal{F}(\mathbf{f})}\right\rVert}$,
where ${\left\lVert{\cdot}\right\rVert}$ is the complex Euclidean norm.

###### Proof.

(a) follows directly from the definition. For (b), to show
${\left\lVert{\mathbf{f}}\right\rVert}^{2}={\left\lVert{\mathcal{F}(\mathbf{f})}\right\rVert}^{2}$,
we compute

$\displaystyle{\left\lVert{\mathcal{F}(\mathbf{f})}\right\rVert}^{2}=\sum_{1\leq
j\leq n}\sum_{1\leq i\leq n}|\widehat{\mathbf{f}}_{i}p_{i,j}|^{2}$
$\displaystyle=$ $\displaystyle\sum_{1\leq i\leq
n}|\widehat{\mathbf{f}}_{i}|^{2}\sum_{1\leq j\leq n}|p_{i,j}|^{2}=\sum_{1\leq
i\leq
n}|\widehat{\mathbf{f}}_{i}|^{2}{\left\lVert{\mathbf{v}_{2,i}}\right\rVert}^{2}$
$\displaystyle=$ $\displaystyle\sum_{1\leq i\leq
n}|\widehat{\mathbf{f}}_{i}|^{2}={\left\lVert{\mathbf{f}}\right\rVert}^{2}.$

∎

The proposed framework recovers that of traditional GSP. To see this, suppose
the GSO $\mathbf{S}$ is normal, e.g., symmetric or orthogonal. In this case,
both $\mathbf{V}_{1}$ and $\mathbf{V}_{2}$ can be chosen to be equal. This
implies that $p_{i,j}=\delta_{i,j}=1$ if $i=j$ and $0$ otherwise. Therefore,
$\mathcal{F}(\mathbf{f})_{i,j}=\widehat{\mathbf{f}}_{i}$ if $i=j$ and $0$
otherwise. This is the reason why for undirected graphs or directed cycle
graphs, we only need $n$ (diagonal) components in GSP.

######  Example 1.

In this example, we show sample spectral plots (in Fig. 1) of the same signal
on different directed graphs. Let $G_{0},\ldots,G_{5}$ be directed graphs on
$50$ nodes. For the extreme cases, $G_{0}$ is the undirected cycle graph and
$G_{5}$ is the directed cycle graph. Intermediately, for $0<k<5$, $G_{k}$ is
obtained from $G_{0}$ by making randomly chosen $10k$ edges directed according
to the ordering of the nodes. Shift operators are the (directed) Laplacians
$\mathbf{L}_{k},0\leq k\leq 5$. We choose a bandlimited signal $\mathbf{f}$
w.r.t. $\mathbf{L}_{0}$. The spectral plots of the componentwise absolute
values of its Fourier transforms w.r.t. $\mathbf{L}_{k},0\leq i\leq 5$ are
shown in Fig. 1. For $\mathbf{L}_{0}$, the result is the same as that
traditional GSP offers (placed along the diagonal). If we let $k$ increase, we
see that the patterns change gradually. We shall see that this supports the
result on graph perturbation in Section V below. If the graph differs more
from $G_{0}$, then the spectrum of $\mathbf{f}$ is more spread out, suggesting
the usefulness of having a lossless signal processing framework.

Fig. 1: Spectral plots of $\mathbf{f}$ w.r.t. $\mathbf{L}_{i},0\leq i\leq 5$.

##  IV Convolution

Important tools in GSP are the convolutions [1, 23, 24]. One way to interpret
the convolution is to consider it as a polynomial in the GSO $\mathbf{S}$. We
can do so for GSP for the particular reason that
$\mathbf{S}=\mathbf{W}\bm{\Sigma}\mathbf{W}^{{}^{\intercal}}$ is
diagonalizable and a monomial takes the form
$\mathbf{S}^{m}=\mathbf{W}\bm{\Sigma}^{m}\mathbf{W}^{{}^{\intercal}}$.
Therefore, for the expression $\bm{\Sigma}^{m}$, we essentially take entrywise
multiplication in the frequency domain. However, this does not hold for an
$\mathbf{S}$ from a directed graph, if it is not diagonalizable. Therefore, we
may consider the interpretation that convolution is a multiplication in the
frequency domain [24]. For this, we use the identification of the frequency
domain with $\mathbb{C}^{n^{2}}$. More precisely, we define:

######  Definition 2.

For $\mathbf{g},\mathbf{h}\in\mathbb{C}^{n^{2}}$, let
$\mathbf{g}\odot\mathbf{h}\in\mathbb{C}^{n^{2}}$ be the entrywise
multiplication, i.e., the Hadamard product. Moreover, defined
$\mathbf{h}^{k}=\mathbf{h}\odot\cdots\odot\mathbf{h}$ to be the $k$-fold self
multiplication. For any $\mathbf{h}\in\mathbb{C}^{n^{2}}$, it defines a
_convolution_ $\mathbf{h}*:\mathbb{R}^{n}\to\mathbb{C}^{n}$ by

$\displaystyle\mathbf{h}*\mathbf{f}=\mathcal{I}\Big{(}\widetilde{\mathbf{f}}\odot\mathbf{h}\Big{)}.$

For the polar decomposition of $\mathbf{S}$, let the diagonal entries of
$\bm{\Lambda}_{1}$ (resp. $\bm{\Lambda}_{2}$) be $(\lambda_{1,i})_{1\leq i\leq
n}$ (resp. $(\lambda_{2,i})_{1\leq i\leq n}$).

######  Lemma 2.

Let $\mathbf{h}_{\mathbf{S}}\in\mathbb{C}^{n^{2}}$ be defined by
$\mathbf{h}_{\mathbf{S},i,j}=\lambda_{2,i}\lambda_{1,j}$. Then
$\mathbf{h}_{\mathbf{S}}*\mathbf{f}=\mathbf{S}(\mathbf{f})$.

###### Proof.

From the definition, we directly evaluate $\mathbf{h}_{\mathbf{S}}*\mathbf{f}$
as

$\displaystyle\mathbf{h}_{\mathbf{S}}*\mathbf{f}=\sum_{1\leq j\leq
n}\Big{(}\sum_{1\leq i\leq
n}\widehat{\mathbf{f}}_{i}p_{i,j}\lambda_{2,i}\lambda_{1,j}\Big{)}\mathbf{v}_{1,j}$
$\displaystyle=$ $\displaystyle\sum_{1\leq j\leq
n}\lambda_{1,j}\Big{(}\sum_{1\leq i\leq
n}\widehat{\mathbf{f}}_{i}p_{i,j}\lambda_{2,i}\Big{)}\mathbf{v}_{1,j}$
$\displaystyle=$ $\displaystyle\sum_{1\leq j\leq
n}\lambda_{1,j}\langle\mathbf{v}_{1,j},\sum_{1\leq i\leq
n}\widehat{\mathbf{f}}_{i}\lambda_{2,i}\mathbf{v}_{2,i}\rangle\mathbf{v}_{1,j}$
$\displaystyle=$ $\displaystyle\sum_{1\leq j\leq
n}\lambda_{1,j}\langle\mathbf{v}_{1,j},\mathbf{V}_{2}\bm{\Lambda}_{2}\mathbf{V}_{2}^{H}\mathbf{f}\rangle\mathbf{v}_{1,j}$
$\displaystyle=$ $\displaystyle\sum_{1\leq j\leq
n}\lambda_{1,j}\langle\mathbf{v}_{1,j},\mathbf{P}(\mathbf{f})\rangle\mathbf{v}_{1,j}$
$\displaystyle=$
$\displaystyle\mathbf{U}\big{(}\mathbf{P}(\mathbf{f})\big{)}=\mathbf{S}(\mathbf{f}).$

∎

As $\mathbf{U}$ is unitary and $\mathbf{P}$ is positive semi-definite,
$\lambda_{2,i}$ is always non-negative and it can be interpreted as the
radius, while $\lambda_{1,j}$ is of the form $e^{\mathbf{i}\theta}$. Hence,
$\lambda_{2,i}\lambda_{1,j}$ is the polar form of the complex number
$\mathbf{h}_{i,j}$. As a consequence of Lemma 2, if $\mathbf{S}$ is normal,
then $\mathbf{h}_{\mathbf{S}}^{k}*$ is just the monomial convolution filter
$\mathbf{S}^{k}$.

We have the freedom to choose $\mathbf{h}$ to construct different convolution
filters. For example, if $\mathbf{h}$ is a characteristic signal valued in
$\\{0,1\\}$, then we have a bandpass filter. Sampling theory can be developed
analogously to [8]. For another construction, consider any function
$f:\mathbb{C}\to\mathbb{C}$. We use $f(\mathbf{h}_{\mathbf{S}})$ to denote the
signal in $\mathbb{C}^{n^{2}}$ whose $(i,j)$-th entry is
$f(\mathbf{h}_{\mathbf{S},i,j})=h(\lambda_{2,i}\lambda_{1,j})$. Examples of
$f$ include polynomials or functions of the form $1/(1+cx)$ (cf. Section VII).

##  V Graph perturbation

In this section, we discuss graph perturbation, which is used in [19] to
modify a directed graph to obtain a diagonalizable shift operator. Many
procedures are considered as a perturbation. For example, adding or removing a
small percentage of directed edges is considered as a perturbation. Another
example is a small change in the edge weights if the graph is weighted. Though
our framework is lossless without modifying the graph structure, we still want
to understand how such perturbations affect the Fourier transform and
filtering.

In this section, we use ${\left\lVert{\cdot}\right\rVert}$ to denote the
matrix operator norm, and ${\left\lVert{\cdot}\right\rVert}_{F}$ to denote the
matrix Frobenius norm, the same as the $2$-norm. All norms are equivalent for
finite dimensional spaces. For example, for any matrix $\mathbf{M}$ of rank
$r$, we have
${\left\lVert{\mathbf{M}}\right\rVert}\leq{\left\lVert{\mathbf{M}}\right\rVert}_{F}\leq\sqrt{r}{\left\lVert{\mathbf{M}}\right\rVert}$.

For another shift operator $\mathbf{S}^{\prime}$, denote
$\mathbf{S}^{\prime}-\mathbf{S}$ by $\delta\mathbf{S}^{\prime}$, as the
perturbation matrix. Usually, we are interested in small perturbation, which
means ${\left\lVert{\delta\mathbf{S}^{\prime}}\right\rVert}$ is small. In this
subsection, we assume that $\mathbf{S}$ and $\mathbf{S}^{\prime}$ are
invertible and _non-derogatory_ , i.e., both factors in the polar
decomposition do not have repeated eigenvalues. If
${\left\lVert{\delta\mathbf{S}^{\prime}}\right\rVert}$ is sufficiently small,
then $\mathbf{S}$ being invertible and non-derogatory implies the same holds
for $\mathbf{S}^{\prime}$.

To state the following observation, let $P_{m}(x,y)$ be the polynomial
$(x+y)^{m}-y^{m}$. We notice that each of its monomial terms has at least one
$x$ factor.

######  Lemma 3.

Suppose sequences of matrices $\mathbf{M}_{1},\ldots,\mathbf{M}_{m}$ and
$\mathbf{M}_{1}^{\prime},\ldots,\mathbf{M}_{m}^{\prime}$ satisfy: for $1\leq
i\leq m$, matrices $\mathbf{M}_{i},\mathbf{M}_{i}^{\prime}$ are of the same
size and $\mathbf{M}=\prod_{1\leq i\leq m}\mathbf{M}_{i}$,
$\mathbf{M}^{\prime}=\prod_{1\leq i\leq m}\mathbf{M}_{i}^{\prime}$ are well-
defined. Let
$\delta\mathbf{M}_{i}^{\prime}=\mathbf{M}_{i}^{\prime}-\mathbf{M}_{i}$. If
there are $x,y>0$ such that
${\left\lVert{\delta\mathbf{M}_{i}^{\prime}}\right\rVert}\leq
x,{\left\lVert{\mathbf{M}_{i}}\right\rVert}\leq y,1\leq i\leq m$, then
${\left\lVert{\mathbf{M}^{\prime}-\mathbf{M}}\right\rVert}\leq P_{m}(x,y)$.

###### Proof.

Notice that
${\left\lVert{\mathbf{M}^{\prime}-\mathbf{M}}\right\rVert}={\left\lVert{\prod_{1\leq
i\leq m}(\mathbf{M}_{i}+\delta\mathbf{M}_{i}^{\prime})-\prod_{1\leq i\leq
m}\mathbf{M}_{i}}\right\rVert}$. We expand the product $\prod_{1\leq i\leq
m}(\mathbf{M}_{i}+\delta\mathbf{M}_{i}^{\prime})$ and apply the triangle
inequality. The norm
${\left\lVert{\mathbf{M}^{\prime}-\mathbf{M}}\right\rVert}$ is upper bounded
by a sum of the form ${\left\lVert{\prod_{1\leq i\leq
m}\mathbf{N}_{i}}\right\rVert}$ with $\mathbf{N}_{i}$ either $\mathbf{M}_{i}$
or $\delta\mathbf{M}_{i}^{\prime}$, and at least one of them is
$\delta\mathbf{M}_{i}^{\prime}$. Then as ${\left\lVert{\cdot}\right\rVert}$ is
submultiplicative, we obtain the desired an upper bound $P_{m}(x)$ by
replacing $\mathbf{M}_{i}$ with $y$ and $\delta\mathbf{M}_{i}^{\prime}$ by
$x$. ∎

The Fourier transform Definition 1 and convolutions are linear
transformations, and hence it makes sense to talk about their operator norms
and Frobenius norms.

######  Theorem 1.

There are constants $\epsilon$ and $C$ independent of
$\delta\mathbf{S}^{\prime}$ such that if
${\left\lVert{\delta\mathbf{S}^{\prime}}\right\rVert}_{F}\leq\epsilon$, then
the following holds

  1. (a)

${\left\lVert{\mathcal{F}_{\mathbf{S}}-\mathcal{F}_{\mathbf{S}^{\prime}}}\right\rVert}\leq
P_{3}(C{\left\lVert{\delta\mathbf{S}^{\prime}}\right\rVert}_{F},1)$.

  2. (b)

${\left\lVert{\mathbf{h}_{\mathbf{S}}^{k}*-\mathbf{h}_{\mathbf{S}^{\prime}}^{k}*}\right\rVert}\leq
P_{2k+4}(C{\left\lVert{\delta\mathbf{S}^{\prime}}\right\rVert}_{F},\max(1,\lambda_{2,n}))$.

###### Proof.

Let $\mathbf{V}_{i}^{\prime},\bm{\Lambda}_{i}^{\prime},i=1,2$ be the unitary
and diagonal factors resulting from the polar decomposition of
$\mathbf{S}^{\prime}$. Denote
$\delta\mathbf{V}_{i}^{\prime}=\mathbf{V}_{i}^{\prime}-\mathbf{V}_{i},\delta\bm{\Lambda}_{i}^{\prime}=\bm{\Lambda}_{i}^{\prime}-\bm{\Lambda}_{i},i=1,2$.
We first claim that there is a constant $C$ independent of
$\mathbf{S}^{\prime}$ and $n$ such that for $i=1,2$, then the following holds

  1. (i)

${\left\lVert{\delta\mathbf{V}_{i}^{\prime}}\right\rVert}_{F}\leq
C{\left\lVert{\delta\mathbf{S}^{\prime}}\right\rVert}_{F}$, and

  2. (ii)

${\left\lVert{\delta\bm{\Lambda}_{i}^{\prime}}\right\rVert}_{F}\leq
C{\left\lVert{\delta\mathbf{S}^{\prime}}\right\rVert}_{F}$.

The inequalities follow from the matrix perturbation theory as follows. Let
$\mathbf{S}^{\prime}=\mathbf{U}^{\prime}\mathbf{P}^{\prime}$ be the unique
polar decomposition of the invertible matrix $\mathbf{S}^{\prime}$ and write
$\delta\mathbf{U}^{\prime}=\mathbf{U}^{\prime}-\mathbf{U},\delta\mathbf{P}^{\prime}=\mathbf{P}^{\prime}-\mathbf{P}$.
Then by [25, Theorem 1] and [26], we have
${\left\lVert{\delta\mathbf{U}^{\prime}}\right\rVert}_{F}\leq
c_{0}{\left\lVert{\delta\mathbf{S}^{\prime}}\right\rVert}_{F}$ and
${\left\lVert{\delta\mathbf{P}^{\prime}}\right\rVert}_{F}\leq\sqrt{2}{\left\lVert{\delta\mathbf{S}^{\prime}}\right\rVert}_{F}$,
where $c_{0}$ depends only on the singular values of $\mathbf{S}$.

By Mirsky’s theorem ([27, Theorem 2]), we have
${\left\lVert{\delta\bm{\Lambda}_{2}^{\prime}}\right\rVert}_{F}\leq{\left\lVert{\delta\mathbf{P}^{\prime}}\right\rVert}_{F}$.
On the other hand, by the Hoffman-Wielandt theorem ([28, VI.34]) for normal
matrices, there is a permutation of $\sigma$ of indices $\\{1,\ldots,n\\}$
such that

$\displaystyle\sum_{1\leq i\leq
n}|\lambda_{1,i}-\lambda_{1,\sigma(i)}^{\prime}|^{2}\leq{\left\lVert{\delta\mathbf{U}^{\prime}}\right\rVert}_{F}^{2}\leq
c_{0}^{2}{\left\lVert{\delta\mathbf{S}^{\prime}}\right\rVert}_{F}^{2}\leq
c_{0}^{2}\epsilon^{2}.$

As all $\lambda_{1,i},1\leq i\leq n$ are distinct by the non-derogatory
condition, if $\epsilon$ is sufficiently small, then the closet eigenvalue to
$\lambda_{1,i}$ is $\lambda^{\prime}_{1,i}$. Hence, $\sigma$ is the identity
permutation and we have
${\left\lVert{\delta\bm{\Lambda}_{1}^{\prime}}\right\rVert}_{F}\leq{\left\lVert{\delta\mathbf{U}^{\prime}}\right\rVert}_{F}\leq
c_{0}\epsilon$. Therefore, together with the assumption that both $\mathbf{U}$
and $\mathbf{P}$ are non-deregatory, we may choose $\epsilon$ small enough
such that there is a $\rho>0$ (depending only on $\epsilon$) such that
$|\lambda_{1,i}-\lambda_{1,j}^{\prime}|\geq\rho,|\lambda_{2,i}-\lambda_{2,j}^{\prime}|\geq\rho$
for any $1\leq i\neq j\leq n$.

By the Davis-Kahan theorem ([28, VII.3]) and its generalization (for normal
matrices) in [29, Corollary 5], there is a constant $c_{1}$ depending only on
$\rho$ such that
${\left\lVert{\delta\mathbf{V}_{1}^{\prime}}\right\rVert}_{F}\leq
c_{1}{\left\lVert{\delta\mathbf{U}^{\prime}}\right\rVert}_{F}$,
${\left\lVert{\delta\mathbf{V}_{2}^{\prime}}\right\rVert}_{F}\leq
c_{1}{\left\lVert{\delta\mathbf{P}^{\prime}}\right\rVert}_{F}$. Setting
$C=\max(c_{0},\sqrt{2})c_{1}$, we obtain the claimed inequalities. From the
proof, we see that $C$ depends only on the singular values of $\mathbf{S}$ and
$\epsilon$, which depends only on the eigenvalues of $\mathbf{U}$ and
$\mathbf{P}$. This proves the claim.

To prove (a) and (b), we apply Lemma 3 using the claim above. For (a), recall
that by Lemma 1, we have
$\mathcal{F}_{\mathbf{S}}(\mathbf{f})=\mathbf{V}_{1}^{H}\mathbf{V}_{2}D(\mathbf{V}_{2}^{H}\mathbf{f})$.
Notice that $D(\mathbf{V}_{2}^{H}\mathbf{f})$ is a linear operator with a
$n^{2}$-dimensional codomain. Its norm is $1$ and satisfies the same
perturbation bound as $\mathbf{V}_{2}^{H}$. Therefore (a) follows immediately
from Lemma 3. For (b), we observe that $\mathbf{h}_{\mathbf{S}}^{k}*$ can be
decomposed into $2k+4$ factors:

  * •

$D(\mathbf{V}_{2}^{H}\mathbf{f})$, $\mathbf{V}_{2}$, $\mathbf{V}_{1}^{H}$ from
$\mathcal{F}_{\mathbf{S}}$;

  * •

$\mathbf{V}_{1}$ from $\mathcal{I}_{\mathbf{S}}$;

  * •

$k$-copies of the transformation $\mathbb{C}^{n^{2}}\to\mathbb{C}^{n^{2}}$
that multiplies the $j$-th column by $\lambda_{1,j},1\leq j\leq n$.

  * •

$k$-copies of the transformation $\mathbb{C}^{n^{2}}\to\mathbb{C}^{n^{2}}$
that multiplies the $i$-th row by $\lambda_{2,i},1\leq i\leq n$.

Each of these factors is bounded in norm by $\max(1,\lambda_{2,n})$ and its
perturbation bound is
$C{\left\lVert{\delta\mathbf{S}^{\prime}}\right\rVert}_{F}$ by the above
claims. Therefore, (b) of the theorem follows from Lemma 3. ∎

Intuitively, the result claims that if we approximate $\mathbf{S}$ by a normal
$\mathbf{S}^{\prime}$ such that traditional GSP tools can be applied, then the
framework developed in this paper “almost” agrees with the traditional GSP.

##  VI Miscellaneous

This section contains a few miscellaneous topics.

###  VI-A Different choices of eigenbasis

The polar decomposition $\mathbf{S}=\mathbf{U}\mathbf{P}$ may not be unique.
Moreover, even for fixed $\mathbf{U}$ or $\mathbf{P}$, they may have repeated
eigenvalues, thus violating the non-derogatory condition of Theorem 1. As a
consequence, there is an ambiguity in the matrices
$\mathbf{V}_{1},\mathbf{V}_{2}$ needed to define the Fourier transform
Definition 1. To (partially) resolve this issue, we propose an optimization
framework.

Let $d(\cdot,\cdot)$ be a metric on the space of $n\times n$ unitary matrices
$\mathcal{U}_{n}$ such as the metric induced by the operator norm. For a polar
decomposition $\mathbf{S}=\mathbf{U}\mathbf{P}$, we propose to construct the
unitary matrices of eigenvectors $\mathbf{V}_{1},\mathbf{V}_{2}$ by solving
the following optimization

$\displaystyle\min_{\mathbf{V}_{1},\mathbf{V}_{2}\in\mathcal{U}_{n}}d(\mathbf{V}_{1},\mathbf{V}_{2}),\text{
provided }\mathbf{V}_{1},\mathbf{V}_{2}\text{ form eigenbasis of
}\mathbf{U},\mathbf{P}.$ (5)

An important example is when $\mathbf{S}$ is the Laplacian of a connected
undirected graph $G$. In this case, by convention, we apply the canonical
polar decomposition $\mathbf{S}=\mathbf{I}\mathbf{S}$, where $\mathbf{I}$ is
the identity matrix. Solving the above optimization results in
$\mathbf{V}_{1}=\mathbf{V}_{2}$ and $d(\mathbf{V}_{1},\mathbf{V}_{2})=0$. This
agrees with what is being considered in GSP. Suppose $\mathbf{S}^{\prime}$ is
the Laplacian of connected undirected $G^{\prime}$ that is a small
perturbation of $G$. Then we still have a constant $C$ independent of
$G^{\prime}$ such that
${\left\lVert{\mathcal{F}_{\mathbf{S}}-\mathcal{F}_{\mathbf{S}^{\prime}}}\right\rVert}\leq
C{\left\lVert{\delta\mathbf{S}^{\prime}}\right\rVert}_{F}$ as long as
$\mathbf{S}$ does not have repeated eigenvalues by Theorem 1. Similarly, if
$\mathbf{S}$ is the adjacency matrix of the directed cycle graph, then for the
canonical polar decomposition $\mathbf{S}=\mathbf{S}\mathbf{I}$, we again have
$\mathbf{V}_{1}=\mathbf{V}_{2}$ by solving (5) with
$d(\mathbf{V}_{1},\mathbf{V}_{2})=0$.

More generally, the same discussion works if $\mathbf{S}$ is _normal_ , in
which case $\mathbf{U},\mathbf{P}$ are simultaneously diagonalizable, and
solving (5) also yields $\mathbf{V}_{1}=\mathbf{V}_{2}$ with
$d(\mathbf{V}_{1},\mathbf{V}_{2})=0$.

###  VI-B An overall picture

For simplicity, our discussion in Section V focuses on the space $\mathcal{S}$
of all invertible and non-derogatory shift operators, which can be weighted.
The space $\mathcal{S}$ is endued naturally with the metric
$d_{\mathcal{S}}(\cdot,\cdot)$ induced by the operator norm
${\left\lVert{\cdot}\right\rVert}$. In this paper, we have defined for each
$\mathbf{S}\in\mathcal{S}$ the Fourier transform operator
$\mathcal{F}_{\mathbf{S}}$. Moreover, for each $k$, we have the convolution
operator $\mathbf{h}_{\mathbf{S}}^{k}*$.

The traditional GSP theory is incomplete by defining GFT (and convolutions) on
the subspace of $\mathcal{S}$ consisting of symmetric shift operators. In this
paper, we have extended the constructions to all of $\mathcal{S}$, which
subsumes GSP as a special case. Moreover, by Theorem 1, the extension of the
traditional theory is in a continuous manner. We call our framework _directed
GSP_ (DGSP). A comparison with GSP is summarized in Table I.

TABLE I: Comparison between traditional GSP and the proposed approach “DGSP” | GSP | DGSP
---|---|---
Applicability | Symmetric matrices | $\mathcal{S}$
Frequency domain | $\mathbb{R}^{n}$ | $\mathbb{C}^{n^{2}}$
Fourier transform | $\mathbf{V}^{{}^{\intercal}}\mathbf{f}$ | $\mathbf{V}_{1}^{H}\mathbf{V}_{2}D(\mathbf{V}_{2}^{H}\mathbf{f})$
Inverse transform | $\mathbf{V}\mathbf{f}$ | $\mathbf{V}_{1}\mathbf{M}\mathbf{c}$
Convolution | Polynomial in $\mathbf{S}$ | Not polynomial in $\mathbf{S}$
Relation to GSP | $--$ | GSP if $\mathbf{S}$ is normal

###  VI-C General matrix decomposition

The Fourier transform (4) has the element of the chain rule in calculus. More
precisely, for a signal $\mathbf{f}$, the polar decomposition trivially
implies that $\mathbf{S}\mathbf{f}=\mathbf{U}(\mathbf{P}\mathbf{f})$. In the
expression
$\langle\mathbf{v}_{2,i},\mathbf{f}\rangle\langle\mathbf{v}_{1,j},\mathbf{v}_{2,i}\rangle$,
we view $\langle\mathbf{v}_{2,i},\mathbf{f}\rangle$ as the “partial derivative
of $\mathbf{f}$ w.r.t. $\mathbf{v}_{2,i}$, and the factor
$\langle\mathbf{v}_{1,j},\mathbf{v}_{2,i}\rangle$ as the “partial derivative”
of $\mathbf{v}_{2,i}$ w.r.t. $\mathbf{v}_{1,j}$. This gives us hints on how to
generalize the framework of the paper.

Instead of considering the polar decomposition
$\mathbf{S}=\mathbf{U}\mathbf{P}$, let

$\displaystyle\mathbf{S}=\mathbf{M}_{1}\cdots\mathbf{M}_{k},k>1$ (6)

be any decomposition of $\mathbf{S}$ into a product of normal matrices
$\mathbf{M}_{i},1\leq i\leq k$. Assume that for each $\mathbf{M}_{i}$, there
is an eigendecomposition
$\mathbf{M}_{i}=\mathbf{V}_{i}\bm{\Lambda}_{i}\mathbf{V}_{1}^{H}$ such that
$\mathbf{V}_{i}$ is unitary and $\bm{\Lambda}_{i}$ is the diagonal matrix of
eigenvalues. Let the columns of $\mathbf{V}_{j}$ be $\\{\mathbf{v}_{i,j},1\leq
j\leq n\\}$. Then the Fourier transform of a signal $\mathbf{f}$ is a map
$\mathbb{R}^{n}\to\mathbb{C}^{n^{k}}$ defined by

$\displaystyle\mathbf{f}\mapsto\Big{(}\prod_{1\leq i\leq
k-1}\langle\mathbf{v}_{i,j_{i}},\mathbf{v}_{i+1,j_{i+1}}\rangle\Big{)}\langle\mathbf{v}_{k,j_{k}},\mathbf{f}\rangle,1\leq
j_{1},\ldots,j_{k}\leq n,$

as a generalization of (4) in Definition 1.

If the eigendecompositions of $\mathbf{M}_{i}$’s are non-unique, as in the
case of the polar decomposition, we may jointly estimate $\mathbf{V}_{i},1\leq
i\leq k$ by solving the following optimization taking the form analogous to
(5):

$\displaystyle\min_{\mathbf{V}_{1},\ldots,\mathbf{V}_{k}}\sum_{1\leq i\leq
k-1}d(\mathbf{V}_{i},\mathbf{V}_{i+1}),\text{ provided }\mathbf{V}_{i}\text{
forms an eigenbasis of }\mathbf{M}_{i}.$

We notice that the codomain of the Fourier transform is $n^{k}$ dimensional.
Therefore, it will be intractable if $k$ is large. This means that it is only
reasonable to consider only the matrix decomposition of a small number of
factors. Examples of such a matrix decomposition include the reverse polar
decomposition, the singular value decomposition (SVD) [22], and the Sinkhorn
decomposition involving doubly stochastic matrix [30].

##  VII Numerical results

In this section, we consider the heat flow dataset of the Intel Berkeley
Research lab.111http://db.csail.mit.edu/labdata/labdata.html Temperature data
are collected from $53$ sensors, denoted by $V$, placed in the lab. We use the
setup of [18] and construct a planar graph $G$ with an edge set $E$ of size
$87$. Each edge is given a direction such that heat flows from high
temperature to low temperature. A signal $\mathbf{f}$ on $G$ consists of
temperature readings on the sensors $V$.

We follow [18, Section 4.2] and perform the signal denoising task. For a
temperature signal $\mathbf{f}$, we add a noise $\mathbf{n}$ whose components
are independent Gaussians with mean $\mu=0$ and standard deviation $\sigma$.
The task is to recover $\mathbf{f}$, while observing $\mathbf{f}+\mathbf{n}$.
In Fig. 2, we show the spectral plots of a true signal $\mathbf{f}$ and a
noise $\mathbf{n}$. We see that the spectral of $\mathbf{f}$ is mostly
supported for small frequency indices as expected. On the other hand, high
magnitude components of the spectral coefficients $\mathbf{n}$ are spread over
the entire frequency domain. This suggests a low-pass filter can be used for
denoising. As in [18], we consider a convolution filter
$f(\mathbf{h}_{\mathbf{S}})$ (Section IV) with $f(x)=1/(1+cx)$ for a tunable
hyperparameter $c$.

Fig. 2: The sample spectral plots of a temperature signal in the lab and a
noise signal with $\sigma=4$. The magnitudes of the spectral coefficients are
shown.

Let $\mathbf{f}^{\prime}$ be the recovered signal. We evaluate the performance
by computing the _root-mean-squared-error_ (RMSE)
${\left\lVert{Re(\mathbf{f}^{\prime})-\mathbf{f}}\right\rVert}_{2}/\sqrt{53}$
over all the sensors. Here $Re(\mathbf{f}^{\prime})$ is the componentwise real
part of $\mathbf{f}^{\prime}$. We call our approach “DGSP”. We compare with
the method “Hermitian” of [18]. The boxplots of RMSEs over $1000$ runs and
$\sigma=1,\ldots,8$ are shown in Fig. 3. For reference, we also show the RMSE
between $\mathbf{f}+\mathbf{n}$ and $\mathbf{f}$, called “Base”. We see that
our lossless directed graph signal processing approach “DGSP” performs better
for all $\sigma$.

Fig. 3: The RMSE of noisy signal recovery.

##  VIII Conclusions

This paper introduces a novel signal processing framework for directed graphs.
Our approach utilizes the polar decomposition to introduce a novel Fourier
transform with a larger frequency domain than that of traditional GSP. This
enables us to define convolution using a standard procedure. Our method
possesses two notable features: it ensures losslessness, as the shift operator
can be entirely reconstructed from the polar decomposition factors, and it
encompasses traditional graph signal processing when applied to directed
graphs. For future work, we shall study how the framework might be used for
learning directed graph topology.

## References

  * [1] D. Shuman, S. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” _IEEE Signal Process. Mag._ , vol. 30, no. 3, pp. 83–98, 2013.
  * [2] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on graphs,” _IEEE Trans. Signal Process._ , vol. 61, no. 7, pp. 1644–1656, 2013.
  * [3] M. Tsitsvero, S. Barbarossa, and P. Di Lorenzo, “Signals on graphs: Uncertainty principle and sampling,” _IEEE Trans. Signal Process._ , vol. 64, no. 18, pp. 4845–4860, 2016.
  * [4] A. Ortega, P. Frossard, J. Kovacevic, J. Moura, and P. Vandergheynst, “Graph signal processing: Overview, challenges, and applications,” _Proc. IEEE_ , vol. 106, no. 5, pp. 808–828, 2018.
  * [5] F. Ji and W. P. Tay, “A Hilbert space theory of generalized graph signal processing,” _IEEE Trans. Signal Process._ , vol. 67, no. 24, pp. 6188 – 6203, 2019.
  * [6] G. Leus, A. Marques, J. Moura, A. Ortega, and D. Shuman, “Graph signal processing: History, development, impact, and outlook,” _IEEE Signal Process. Mag._ , vol. 40, no. 4, pp. 49–60, Jun. 2023.
  * [7] X. Jian, F. Ji, and W. P. Tay, “Generalizing graph signal processing: High dimensional spaces, models and structures,” _Foundations and Trends® in Signal Processing_ , vol. 17, no. 3, pp. 209–290, 2023.
  * [8] F. Ji, W. P. Tay, and A. Ortega, “Graph signal processing over a probability space of shift operators,” _IEEE Trans. Signal Process._ , vol. 71, pp. 1159–1174, 2023.
  * [9] Y. Sarcheshmehpour, M. Leinonen, and A. Jung, “Federated learning from big data over networks,” in _ICASSP_ , 2021, pp. 3055–3059.
  * [10] J. Giraldo and T. Bouwmans, “Semi-supervised background subtraction of unseen videos: Minimization of the total variation of graph signals,” in _ICIP_. IEEE, 2020, pp. 3224–3228.
  * [11] A. Sandryhaila and J. Moura, “Discrete signal processing on graphs: Frequency analysis,” _IEEE Trans. Signal Process._ , vol. 62, no. 12, pp. 3042–3054, 2014.
  * [12] J. Domingos and J. M. F. Moura, “Graph Fourier transform: A stable approximation,” _IEEE Trans. Signal Process._ , vol. 68, pp. 4422–4437, 2020\.
  * [13] Y. Chen, C. Cheng, and Q. Sun, “Graph Fourier transform based on singular value decomposition of directed laplacian,” 2022.
  * [14] V. Satuluri and S. Parthasarathy, “Symmetrizations for clustering directed graphs,” in _EDBT_ , 2011, p. 343–354.
  * [15] S. Sardellitti and S. Barbarossa, “On the graph Fourier transform for directed graphs,” _IEEE J Sel Top Signal Process._ , vol. 11, no. 6, pp. 796–811, 2017.
  * [16] H. Mhaskar, “A unified framework for harmonic analysis of functions on directed graphs and changing data,” _Appl Comput Harmon Anal._ , vol. 44, no. 3, pp. 611–644, 2018.
  * [17] B. S. Dees, L. Stankovic, M. Dakovic, A. G. Constantinides, and D. P. Mandic, “Unitary shift operators on a graph,” 2019.
  * [18] S. Furutani, T. Shibahara, M. Akiyama, K. Hato, and M. Aida, “Graph signal processing for directed graphs based on the Hermitian laplacian,” in _ECML/PKDD_ , 2019.
  * [19] B. Seifert and M. Püschel, “Digraph signal processing with generalized boundary conditions,” _IEEE Trans. Signal Process._ , vol. 69, pp. 1422–1437, 2021.
  * [20] B. Hall, _Lie Groups, Lie Algebras, and Representations: An Elementary Introduction_ , 2nd ed. Springer, 2015.
  * [21] S. Kwak, L. Shimabukuro, and A. Ortega, “Frequency analysis and filter design for directed graphs with polar decomposition,” _arXiv preprint arXiv:2312.11421_ , 2023.
  * [22] T. Tao, _Topics in random matrix theory_. American Mathematical Soc., 2012.
  * [23] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in _NeurIPS_ , 2016.
  * [24] F. Ji, W. P. Tay, and A. Ortega, “The faces of convolution: from the fourier theory to algebraic signal processing,” _arXiv preprint arXiv:2307.07926_ , 2023.
  * [25] R. Li, “New perturbation bounds for the unitary polar factor,” _SIAM J. Matrix Anal. Appl._ , vol. 16, pp. 327–332, 1995.
  * [26] C. Chen and J. Sun, “Perturbation bounds for the polar factors,” _J. Comput. Math._ , vol. 7, no. 4, pp. 397–401, 1989.
  * [27] G. Stewart, “Perturbation theory for the singular value decomposition,” 1990.
  * [28] R. Bhatia, _Matrix Analysis_. Springer Science and Business Media, 1996.
  * [29] S. Bhattacharya, “On some bounds on the perturbation of invariant subspaces of normal matrices with application to a graph connection problem,” _J. Inequal. Appl._ , no. 1, p. 75, 2022.
  * [30] R. Sinkhorn, “A relationship between arbitrary positive matrices and doubly stochastic matrices,” _Ann. Math. Statist._ , vol. 35, p. 876–879, 1964\.


An approach to cold dark matter deviation and the H_0 tension problem
by using machine learning

Emilio Elizalde^1,2[E-mail<EMAIL_ADDRESS>Janusz Gluza^3[E-mail<EMAIL_ADDRESS>Martiros Khurshudyan^1,2,3[E-mail<EMAIL_ADDRESS>
$^{1}$ Consejo Superior de Investigaciones Científicas, ICE/CSIC-IEEC,
Campus UAB, Carrer de Can Magrans s/n, 08193 Bellaterra (Barcelona) Spain

$^{2}$ International Laboratory for Theoretical Cosmology, Tomsk State University of Control Systems and Radioelectronics (TUSUR), 634050 Tomsk, Russia

$^{3}$ Institute of Physics, University of Silesia,  Katowice, Poland

In this work, two different models, one with cosmological constant $\Lambda$, and baryonic and dark matter (with $\omega_{dm} \neq 0$), and the other with an $X$ dark energy (with $\omega_{de} \neq -1$), and baryonic and dark matter (with $\omega_{dm} \neq 0$), are investigated and compared. Using  Bayesian machine learning analysis, constraints on the free parameters of both models are obtained for the three redshift ranges: $z\in [0,2]$, $z\in [0,2.5]$, and $z\in [0,5]$, respectively.
For the first two redshift ranges, high-quality observations of the expansion rate $H(z)$ exist already, and they are used for validating the fitting results. Additionally, the extended range $z\in [0,5]$ provides predictions of the model parameters, verified when reliable higher-redshift $H(z)$ data are available. This learning procedure, based on the expansion rate data generated from the background dynamics of each model, shows that, at cosmological scales, there is a deviation from the cold dark matter paradigm, $\omega_{dm} \neq 0$, for all three redshift ranges. The results show that this approach may qualify as a solution to the $H_{0}$ tension problem. Indeed, it hints at how this issue could be effectively solved (or at least alleviated) in cosmological models with interacting dark energy.

§ INTRODUCTION

Recently, an interesting issue which goes under the name of the $H_{0}$ tension problem has appeared. Specifically, the Planck-CMB data analysis provides a value $H_{0} = 67.4 \pm 0.5$ km/s/Mpc when the $\Lambda$CDM scenario is assumed <cit.>.
But, a local measurement coming from the Hubble Space Telescope yields $H_{0} = 73.52 \pm 1.62$ km/s/Mpc <cit.>.
Thus, it has been pointed out that there is an actual tension between the calculated values of $H_{0}$ obtained from these two observational data sets, see also  <cit.> and <cit.>.
The difference between the two observed values of  $H_{0}$ is unexpectedly huge. Thus, on top of the other problems of modern cosmology, we have an additional one, which indicates that either there is new physics here or maybe a serious problem with either of the two measurements. Indeed, at this moment, it is very hard to estimate where the problem lies since we still need to reveal the nature of dark energy from one side, and the nature of dark matter, from the other. Anyhow, if the problem would eventually point towards new physics, then we would be in a perfect starting situation to clarify, constrain and select good dark energy models as extensions and modifications of General Relativity. However, if it turned out to be a measurement-related issue, then we would be dealing with an artificially crafted problem, not carrying useful information for cosmology. In any case, it is too early to reach a final conclusion, as has been discussed in the recent literature. It has been shown that some of the interacting dynamical dark energy models available on the market could be employed to solve the issue. But also, a modification of General Relativity may well lead to other solutions

(to mention some references).

The idea of considering interacting dark energy models is not new and has been used successfully in cosmology over many years. In interacting dark energy scenarios, one operates at the macroscopic level and introduces a non-gravitational interaction between dark energy and dark matter. This can be viewed as an artificial intervention and, to our knowledge, up to this moment, no well-defined theory has been found that can explain where the non-gravitational interaction comes from. One draws the existence of non-gravitational interactions between dark energy and dark matter merely from observations. However, we should recognize that this could be an artifact created  by a non-proper screening of our observational Universe. This could induce an improper bias, which added to the lack of a full understanding of the accelerated expansion of the Universe, could create the illusion of the additional interaction. However, we would like to stress that, in no way the above discussion wants to diminish the real value of the huge amount of work done in this direction, while exploring such possibility. It is just a motivation for additional research following a different path, since we do believe that there is yet job to be done in order to answer such questions. The community is still very active in exploring new and interesting aspects of interacting dark energy models <cit.>
(to mention some references).

Before focusing on the main goal of our work, we would like to recall a way to infer cosmological knowledge from observational data. It is based on the well known Bayes theorem, which allows connecting the posterior distribution with the prior one with a likelihood. It has been applied very often to constrain the models yielding new knowledge. We shall omit technical details, but we would like to mention that we need to use a very well-defined sampling technique to build the posterior distribution. In the literature, there are interesting approaches mainly based on different ideas on how to deal with rejected and accepted samplings. The last ones significantly improve the statistics and speed up the fitting process. However, we need to have very well-defined models to proceed with the calculations efficiently.

Retaking the above discussion on interacting dark energy models, let us now clarify what we mean by a well-defined model; namely, in the process, we need to supply an explicit form of dark energy, dark matter and the form of the non-gravitational interaction between them. Then, at the end of the day, constraining the model, we can make some inference. In particular, we can claim whether a specific type of non-gravitational interacting dark energy model is supported or not. However, there is another method known as Gaussian Process, which has been used in the recent literature very intensively (see for instance <cit.>
and <cit.>
and references therein for more details). It is one of the Machine Learning approaches allowing to do reconstruction directly from data. In this case, for the same interacting dark energy models, we can reconstruct the possible non-gravitational interaction term as soon as the models of dark energy and dark matter are given. In this case, we do not need to supply any parametric form describing the interaction.

In other words, following this approach, we can get the most general information about the specific task in a model independent way just having access to the observational data. The method has already gained a lot of attention, and we refer the reader to our reference list for more details on the topic. Huge efforts are being made attempting to integrate different Machine Learning algorithms into cosmology, and astrophysics <cit.>.
There are different reasons for this, but the main one is the ability of different Machine Learning algorithms to learn the reason leading to the answer. One of them is the Bayesian (Probabilistic) Machine Learning approach, which, generally speaking, uses probabilities to think about probabilities. Again omitting all details and the history behind the approach, we would like to mention that it is beneficial to gain information about the model, which would be extremely hard to do with the usual approaches requiring data collection. This is because it is based on the generative process directly connected to the model, and, having it, we can explore the regions from which, in practice, it is tough to get answers. It should be mentioned that, recently, this method has been applied with success to study the $H_{0}$ tension problem in a single inhomogeneous viscous Universe. Moreover, it has been demonstrated that the $H_{0}$ tension problem solution is very weakly dependent on the specific form of the viscosity term. Instead, it has been shown that the solution mostly depends on the equation of state parameter of the barotropic fluid, used to craft the viscous dark energy fluid model <cit.>.
On the other hand, using  Bayesian (Probabilistic) Machine Learning, the cosmic opacity has been constrained, revealing its connection with the $H_{0}$ tension problem <cit.>.
In another paper, we surprisingly found that the Swampland criteria itself can be used to constrain the models for the dark energy-dominated Universe <cit.>.
Motivated by the previously obtained results, we here continue to apply the same method to understand the extent to which we can use it and the potential limitations that one could face.

This work aims to reveal possible connections between dark matter and the $H_{0}$ tension problem by using a Bayesian (Probabilistic) Machine Learning approach. The motivation comes from the ability of interacting dark energy models to solve the mentioned tension problem and how non-gravitational interactions may appear in cosmology. In particular, the presence of a non-gravitational interaction will not preserve the energy densities of each component separately anymore, which would force us to accept that, at the cosmological scales, we do not have purely cold dark matter. In this case, it is reasonable to assume that this deviation from cold dark matter can be a factor in solving (or alleviating at least) the $H_{0}$ tension problem. Motivated by the mentioned results, we would like to apply the Bayesian (Probabilistic) Machine Learning approach to study the $H_{0}$ tension problem and see how it is connected to the dark matter issue. It is easy to see that, with interacting dark energy models, the interaction affects the equation of state parameter of both components. It is then tough to separate the dark energy and dark matter contribution in the solution of the $H_{0}$ tension problem.

The models we shall here study are one with cosmological constant $\Lambda$, and baryonic and dark matter (with $\omega_{dm} \neq 0$), and another with an $X$ dark energy (with $\omega_{de} \neq -1$), and baryonic and dark matter (with $\omega_{dm} \neq 0$), which are investigated and compared. Here, $\omega_{dm}$ is considered to be one of the free parameters. The constraints  are obtained for the three redshift ranges: $z\in [0,2]$, $z\in [0,2.5]$, and $z\in [0,5]$, respectively.
The first two redshift ranges are chosen to validate our fit results since we have very good quality observations of the expansion rate $H(z)$ in these cases. Additionally, the extended range $z\in [0,5]$ provides predictions of the model parameters, verified when reliable higher-redshift $H(z)$ data would become available. The learning method based on the generated $H(z)$ expansion rate data from the background dynamics of each model shows that, on cosmological scales, there is a deviation from the cold dark matter paradigm (it is not related to non-gravitational interaction). Moreover, we will encounter a hint that it could solve the tension problem and another on why the parameters of the interaction term should be small, as is discussed in the literature. Of course, to understand to which extend the observed deviation can be used to overcome the various cosmological problems, instead of involving interacting dark energy models, is still the issue, which will have to be tackled down in forthcoming papers.

This paper is organized as follows. The philosophy behind the method used here is discussed in Sect. <ref>. In the same section, we present the models to be studied. The results are discussed in  Sect. <ref>, which is followed by an analysis of the consequences, both interesting for the $H_{0}$ tension and for the accelerated expanding Universe issues. The model-based generative process employed in the Bayesian Learning approach is built from Eqs. (<ref>) and (<ref>), by assuming that $\omega_{dm}$ is one of the free parameters of the model. The model-based generative process and the analysis leading to our final results are performed by using PyMC3 <cit.>. Final conclusions of our analysis are given in Sect. <ref>.

§ THE METHOD AND MODELS

Our analysis is based on the Bayesian (Probabilistic) Machine Learning approach, and in this section, we will briefly describe its key points. It is an interesting method that has been used to tackle various problems. Its advantage is due to the generative process being used, which allows exploring regions that are usually unreachable for some reason. Moreover, recent studies show that using it with reinforcement and inverse reinforcement can significantly speed up the learning process and, at the same time, reduce the costs of collecting data. On the other hand, it is clear that to build the generative process, we need to have a very good mathematical model describing the interesting processes for us. Because now only the generated data will be the key ingredient involved in this learning process. Fortunately, over the years, we have learned how to craft good models interesting for cosmology and astrophysics. Therefore, it is reasonable to apply the Bayesian (Probabilistic) Machine Learning approach to different problems from those fields and make some inference.

In particular, it can be used for those problems, which cannot be fully explored by other means due to the low quality of data. Moreover, since we use a generative process, it is possible to do some forecasting, too, which can help design new missions and observations. In the recent literature, the basic ideas behind Bayesian Machine Learning have been well documented. Moreover, various frameworks allow the organization of Bayesian (Probabilistic) Machine Learning significantly, reducing the coding time, which allows concentrating on the physics of the process. One of such frameworks is the python-based PyMC3 framework, which can be easily integrated with other python-based frameworks, too <cit.>
[We refer the readers to the nice documentation supplied with PyMC3, for more details on the Bayesian (Probabilistic) Machine Learning approach. Moreover, an existing set of interesting tasks and its solutions can generate intuition when and how Bayesian (Probabilistic) Machine Learning can be used.]. It should be noted that it uses Theano <cit.>. The last one is a deep learning python-based library to construct probability distributions and implement cutting edge inference algorithms.

Now, if we want to use Bayesian (Probabilistic) Machine Learning, we need to go through the following three steps:

* We need to define the model to be used to provide a so-called generative process. In our case, it will be the cosmological model.

* We should envisage the data to be the data obtained from the generative process. In our case, it will be the expansion rate data generated from the cosmological model, which we are going to study[We direct the readers to

for more details about the possible difference between Machine Learning and Bayesian Machine Learning that rose namely in this step.].

* We then run the learning algorithm to get a brand new distribution over the model parameters and update our prior belief. In other words, we need the prior distributions over the model's free parameters. This will correspond to our initial beliefs. Moreover, we need to update them each time to get the posteriors. In our case, we use PyMC3 and in-built functions to do all required steps.

Now let us discuss the models that we will study in our paper. In particular, we will consider two cosmological models having specific expansion rate history, obtained from a more general one presented below

\begin{equation}\label{eq:HXCDM}
H(z) = H_{0} E(z,\Omega_{dm}, \Omega_{b}, \omega_{dm}, \omega_{de}),
\end{equation}

\begin{equation}\label{eq:EXCDM}
E(z,\Omega_{dm}, \Omega_{b}, \omega_{dm}, \omega_{de}) = \left [ \Omega_{b} (1+z)^{3} + \Omega_{dm} (1+z)^{3(1+\omega_{dm})}  +  (1- \Omega_{dm} - \Omega_{b}) (1+z)^{3(1+\omega_{de})}\right]^{1/2},
\end{equation}

where $H_{0}$, $\Omega_{b}$ and $\Omega_{dm}$ are the Hubble parameter, and the fraction of the baryonic matter and dark matter, at $z=0$, respectively. On the other hand, $\omega_{dm}$ and $\omega_{de}$ are the equation of state parameters describing dark matter and dark energy, respectively. Moreover, if we consider $\omega_{de} = -1$ and $\omega_{dm} = 0$, we reduce the expansion rate given by Eqs. (<ref>) and (<ref>) to the one corresponding to the standard $\Lambda$CDM model. While, if we consider $\omega_{de} \neq -1$, $\omega_{dm} = 0$ and combine dark matter with baryonic matter, we get the known XCDM model.

Our goal is to constrain two cosmological models using the Bayesian (Probabilistic) Machine Learning approach, where the expansion rate data will be generated. The models are, namely, one with a cosmological constant $\Lambda$, and baryonic and dark matter (with $\omega_{dm} \neq 0$), and another with an $X$ dark energy (with $\omega_{de} \neq -1$), and baryonic and dark matter (with $\omega_{dm} \neq 0$). In both cases, $\omega_{dm}$ is a free parameter and will be constrained in three different redshift ranges, namely $z\in [0,2]$, $z\in [0,2.5]$, and $z\in [0,5]$, respectively. This strategy will allow us to estimate possible deviations of $\omega_{dm}$ from zero and explore how this is related to the $H_{0}$ tension problem. The first two redshift ranges are chosen to validate our fitting results since we have very good quality observations of the expansion rate $H(z)$ in these cases (see Table <ref>).

$z$     $H(z)$     $\sigma_{H}$     $z$     $H(z)$     $\sigma_{H}$

$0.070$     $69$     $19.6$     $0.4783$     $80.9$     $9$

$0.090$     $69$     $12$     $0.480$     $97$     $62$

$0.120$     $68.6$     $26.2$      $0.593$     $104$     $13$

$0.170$     $83$     $8$     $0.680$     $92$     $8$

$0.179$     $75$     $4$      $0.781$     $105$     $12$

$0.199$     $75$     $5$      $0.875$     $125$     $17$

$0.200$     $72.9$     $29.6$      $0.880$     $90$     $40$

$0.270$     $77$     $14$      $0.900$     $117$     $23$

$0.280$     $88.8$     $36.6$      $1.037$     $154$     $20$

$0.352$     $83$     $14$     $1.300$     $168$     $17$

$0.3802$     $83$     $13.5$      $1.363$     $160$     $33.6$

$0.400$     $95$     $17$     $1.4307$     $177$     $18$

$0.4004$     $77$     $10.2$     $1.530$     $140$     $14$

$0.4247$     $87.1$     $11.1$     $1.750$     $202$     $40$

$0.44497$     $92.8$     $12.9$     $1.965$     $186.5$     $50.4$

$$ & $$

$$ & $$

$$ & $$

$0.24$     $79.69$     $2.65$     $0.60$     $87.9$     $6.1$

$0.35$     $84.4$     $7$      $0.73$     $97.3$     $7.0$

$0.43$     $86.45$     $3.68$      $2.30$     $224$     $8$

$0.44$     $82.6$     $7.8$      $2.34$     $222$     $7$

$0.57$     $92.4$     $4.5$      $2.36$     $226$     $8$

$H(z)$ and its uncertainty $\sigma_{H}$ are in the units of km s$^{-1}$ Mpc$^{-1}$. The upper panel consists of thirty samples deduced from the differential age method. The lower panel corresponds to ten samples obtained from the radial BAO method. The table is according to <cit.>
(see also references therein for details).

Additionally, we consider the  $z\in [0,5]$ extended redshift range for forecasting reasons; the full validation of the results will have to wait for the near future when observations of higher redshift $H(z)$ data become available. As we consider the $H_{0}$ tension problem, then the expansion rate data to validate our results is the most obvious choice. The approach can be used to consider other observations, too. However, we wanted to avoid any bias that can be artificially induced from merging different data sets at this stage.

Before ending this section, we should recall the standard notation of FLRW cosmology (with $8 \pi G = c = 1$). In particular, the metric, in this case, has the form

\begin{equation}
ds^{2} = -dt^{2} + a(t)^{2} \sum_{i =1}^{3} (dx^{i})^{2},
\end{equation}

\begin{equation}\label{eq:F1}
H^{2} = \frac{1}{3}\sum_{k}\rho_{k}.
\end{equation}

Moreover, the energy conservation law for each energy source can be expressed as

\begin{equation}\label{eq:EC}
\dot{\rho}_{k} + 3 H (\rho_{k} + P_{k}) = 0,
\end{equation}

with $P_{k} = \omega_{k} \rho_{k}$, where the index $k$ stands for baryons, dark matter and dark energy, respectively. Due to the mentioned assumptions, we end up on the background dynamics given by Eqs. (<ref>) and (<ref>). It should be noted that Eq. (<ref>) gives the evolution of the energy densities of the non-interacting fluids representing the energy budget of the Universe. However, if we consider for instance,  $Q$, a non-gravitational interaction between dark matter and dark energy only, then Eq. (<ref>) turns out to be the following

\begin{equation}\label{eq:EC_b}
\dot{\rho}_{b} + 3 H \rho_{b} = 0,
\end{equation}

\begin{equation}\label{eq:EC_dm}
\dot{\rho}_{dm} + 3 H (\rho_{dm} + P_{dm}) = Q,
\end{equation}

\begin{equation}\label{eq:EC_de}
\dot{\rho}_{de} + 3 H (\rho_{de} + P_{de}) = -Q.
\end{equation}

A closer look to Eqs. (<ref>) and (<ref>) indicates that non-gravitational interaction affects the equation of state parameters of both components. Even if we choose pure cold dark matter, $P_{dm} = 0$, the interaction $Q$ can make it colder or warmer, depending on its form and the constraints on its free parameters. The very brief discussion on interacting dark energy models given above will certainly help the reader follow our results. More details on interacting dark energy models can be found in the references of this paper.

§ RESULTS AND DISCUSSION

In this section, we present and discuss our results. We have divided it into two subsections, intending to simplify the presentation. It had been mentioned already that we consider three redshift ranges for both cosmological models, namely $z\in [0,2]$, $z\in [0,2.5]$ and $z \in [0,5]$, to understand better how our knowledge about the dark matter equation of state parameter can be changed. We remind the reader that our analysis is based on a generative process, where the generated data is the expansion rate of the Universe. This follows from the fact that we are primarily interested in exploring the connection between the $H_{0}$ tension problem and the nature of dark matter. In our opinion, this setup allows us to avoid possible artificial biases that could otherwise appear by using several datasets together, masking the true nature of the problem. Eventually, we can extend this analysis to include other types of "observational" data allowing us to understand the sources of the bias responsible for non-proper conclusions. We will undertake this specific analysis in another paper. The two tables, Table <ref> and Table <ref>, summarize our results here.

§.§ Model with $\omega_{dm} \neq 0$ dark matter and cosmological constant

The first example we have considered is a model of the Universe where the expansion rate is given by Eq. (<ref>) with

\begin{equation}\label{eq:EXCDM_M1}
E(z,\Omega_{dm}, \Omega_{b}, \omega_{dm}) = \left [ \Omega_{b} (1+z)^{3} + \Omega_{dm} (1+z)^{3(1+\omega_{dm})}  +  \Omega_{de} \right]^{1/2},
\end{equation}

where $\Omega_{de} = 1- \Omega_{dm} - \Omega_{b}$. It is a model of the Universe with the cosmological constant as its dark energy and baryonic and dark matter, with $\omega_{dm} \neq 0$. We now devise a generative process based on Bayesian (Probabilistic) Machine Learning to constrain $\omega_{dm}$ and estimate or predict (depending on the redshift range considered) the dynamics of changes that it can undergo based on the high redshift expansion rate measurements. The other parameters of the model to be fitted are $H_{0}$, $\Omega_{dm}$ and $\Omega_{b}$, respectively. On the other hand, as we discuss below, in this way, we can explore other reasons behind the $H_{0}$ tension problem. In particular, someones not related, say, to the dark energy problem. The constraints discussed below result from the learning process based on Eqs. (<ref>) and  (<ref>) using $H_{0} \in [64.0, 80.0]$,  $\Omega_{dm} \in [0.23, 0.4]$, $\Omega_{b} \in [0.0, 0.1]$ and $\omega_{dm} \in [-0.1, 0.1]$ as flat priors. In other words, we have generated data, using Eqs. (<ref>) and (<ref>), by imposing the mentioned priors. The results of the learning procedure can be read from Table <ref> and the ensuing presentation.

* The best fit values of the model parameters are $H_{0} = 68.42 \pm 0.24$ km/s/Mpc, $ \Omega_{dm} = 0.289 \pm 0.007$, $\Omega_{b} = 0.0037 ^{+0.0032}_{+ 0.0025}$, and $\omega_{dm} = -0.035 \pm 0.005$, when $z \in [0,2]$. The contour map is given in Fig. (<ref>), in purple colour.

* The best fit values of the model parameters will most likely be $H_{0} = 68.54 \pm 0.24$ km/s/Mpc,  $ \Omega_{dm} = 0.291 \pm 0.007$, $\Omega_{b} = 0.0035 ^{+ 0.0036}_{-0.0024}$, and $\omega_{dm} = -0.022 \pm 0.005$, when $z \in [0,2.5]$. The contour map is given in Fig. (<ref>), in orange colour.

* The best fit values of the model parameters will most likely be $H_{0} = 68.95 \pm 0.25$ km/s/Mpc,  $ \Omega_{dm} = 0.271 \pm 0.003$, $\Omega_{b} = 0.00035^{+ 0.00038}_{-0.00022}$, and $\omega_{dm} = -0.028 \pm 0.001$, when $z \in [0,5]$. The contour map, when $z \in [0, 5]$, is given in Fig. (<ref>), in red colour.

First, we have investigated the model for $z\in[0,2]$ and found that there is a clear deviation from $\omega_{dm} = 0$. In particular, we have found that its mean value is $\omega_{dm} =  -0.035$. Clearly, the obtained mean for $H_{0}$ is closer to the Planck Mission results than the zero value. On top of this, this slight hints at the possibility that the $H_{0}$ tension problem may be alleviated. On the other hand, the analysis shows that the value of $\omega_{dm}$ will increase when $z\in [0.2.5]$. But still, we have a clear deviation from cold dark matter. Moreover, we have learned that the mean of the free parameters $H_{0}$ and $\Omega_{dm}$  will increase, while the mean value of $\Omega_{b}$ will slightly decrease. The comparison of the results obtained hints here that, most likely, this model cannot fully solve the $H_{0}$ tension problem.  How could it change if we were able to measure the expansion rate up to redshift $z\ = 5$? What would be the result then?

The answer can be seen from Table <ref> which indicates that the mean of the $H_{0}$ parameter will increase compared to the results obtained for the two previous redshift ranges. Moreover, the values of the parameters $\Omega_{dm}$ and $\Omega_{b}$ will decrease dramatically. We also found that $\omega_{dm} \neq 0$ and we have a hint that the deviation from cold dark matter will be preserved. The forecast results associated to the future possible measurements of the expansion rate data of our Universe when $z \in [0,5]$ read as follows: $H_{0} = 68.95 \pm 0.25$ km/s/Mpc,  $\Omega_{dm} = 0.271 \pm 0.003$, $\Omega_{b} = 0.00035^{+ 0.00038}_{-0.00022}$ and $\omega_{dm} = -0.028 \pm 0.001$. To end this subsection, we would like to mention the fact that it seems that the $H_{0}$ tension problem cannot be solved, even if there is a deviation from the usual cold dark matter case. Even if there is a tiny hint supporting the idea, it is actually not enough for the conclusion. Even more, we have observed that the high-redshift expansion rate observations may point to some dynamics of $\omega_{dm}$ and still indicate possible deviations from cold dark matter.

The results discussed in the next subsection show that the assumption taken on the dark energy model can significantly affect the final results. The graphical behaviour of the expansion rate $H(z)$ and the deceleration parameter $q(z)$ corresponding to this second case can be found in Fig. (<ref>). It certifies that we have performed a good learning strategy. Moreover, we see that the learned deviation for $\omega_{dm}$ is not only able to solve the $H_{0}$ tension, but it can also explain the BOSS experiment results for the expansion rate data at $z=2.34$ <cit.>.
This is another interesting outcome of our work here, which requires future and more detailed study. On the other hand, the forecast results for $z\in [0,5]$ allow us to conclude that the high-redshift expansion rate data can validate that the learned $\omega_{dm}$ deviation can additionally be useful to explain the BOSS experiment result. Eventually,  we have learned very tight constraints for the model's free parameters, except for $\Omega_{b}$.

Model 1     $H_{0}$     $\Omega_{dm}$     $\Omega_{b}$     $\omega_{dm}$

when $z\in[0,2]$     $68.42 \pm 0.24$ km/s/Mpc     $0.289 \pm 0.007$     $0.0037 ^{+0.0032}_{+ 0.0025}$      $-0.035 \pm 0.005$

when $z\in[0,2.5]$     $68.54 \pm 0.24$ km/s/Mpc     $0.291 \pm 0.007$     $0.0035 ^{+ 0.0036}_{-0.0024}$      $-0.022 \pm 0.005$

when $z\in[0,5]$      $68.95 \pm 0.25$ km/s/Mpc    $0.271 \pm 0.003$     $0.00035^{+ 0.00038}_{-0.00022}$      $-0.028 \pm 0.001$

Best fit values and $1\sigma$ errors estimated for the $\Lambda$ + Baryonic Matter + Dark Matter (with $\omega_{dm} \neq 0$) model, when $z \in [0,2]$, $z \in [0,2.5]$ and $z \in [0,5]$, respectively. The results have been obtained from a Bayesian (Probabilistic) Machine Learning approach, where the generative process is based on Eq. (<ref>) and Eq. (<ref>) using $H_{0} \in [64.0, 80.0]$,  $\Omega_{dm} \in [0.23, 0.4]$, $\Omega_{b} \in [0.0, 0.1]$ and $\omega_{dm} \in [-0.1, 0.1]$ flat priors. The analysis is based on 10 chains and in each chain, 10,000 "observational" data-sets from the model have been simulated or generated.

\begin{array}{cccc}
\includegraphics[width=120 mm]{model_LCDM.pdf}
\end{array}$

Contour maps of the model given by Eqs. (<ref>) and (<ref>) for $z \in [0,2]$, $z \in [0,2.5]$ and $z \in [0,5]$, respectively. The best fit values of the model parameters have been found to be $H_{0} = 68.42 \pm 0.24$ km/s/Mpc, $ \Omega_{dm} = 0.289 \pm 0.007$, $\Omega_{b} = 0.0037 ^{+0.0032}_{+ 0.0025}$ and $\omega_{dm} = -0.035 \pm 0.005$ when $z \in [0,2]$. While, when  $z \in [0,2.5]$, the best fit values of the model parameters have been found to be $H_{0} = 68.54 \pm 0.24$ km/s/Mpc,  $ \Omega_{dm} = 0.291 \pm 0.007$, $\Omega_{b} = 0.0035 ^{+ 0.0036}_{-0.0024}$ and $\omega_{dm} = -0.022 \pm 0.005$. Eventually, the best fit values of the model parameters have been found to be $H_{0} = 68.95 \pm 0.25$ km/s/Mpc,  $ \Omega_{dm} = 0.271 \pm 0.003$, $\Omega_{b} = 0.00035^{+ 0.00038}_{-0.00022}$ and $\omega_{dm} = -0.028 \pm 0.001$, when $z \in [0,5]$. In all three cases, $H_{0} \in [64.0, 80.0]$,  $\Omega_{dm} \in [0.23, 0.4]$, $\Omega_{b} \in [0.0, 0.1]$, and $\omega_{dm} \in [-0.1, 0.1]$ flat priors have been imposed during the generative process used to generate the "observational" data. The analysis is based on 10 chains and, in each chain, 10,000 "observational" data-sets from the model have been simulated or generated.

\begin{array}{cccc}
\includegraphics[width=80 mm]{H_Model_1.jpeg}&&
\includegraphics[width=80 mm]{q_Model_1.jpeg}\\
\end{array}$

Graphic behaviour of the Hubble parameter compared with known $H(z)$ data is presented on the left-hand side. The purple curve is a plot of the Hubble parameter for the best fit values of the model parameters, when $z \in [0,2]$, the dashed orange curve corresponds to the case when $z \in [0,2.5]$, while the dashed red curve is for the case when $z \in [0,5]$. The red dots represent the known observational $H(z)$ data, which are the same as in Table <ref>. The right-hand side corresponds to the graphic behaviour of the deceleration parameter $q(z)$. The model is given by Eqs. (<ref>) and (<ref>). In both cases, only the best fit values for the model's free parameters obtained by the Bayesian Learning approach and depicted in Table <ref> have been used.

§.§ Model with $\omega_{dm} \neq 0$ dark matter and $\omega_{de} \neq -1$ dark energy

Model 2     $H_{0}$     $\Omega_{dm}$     $\Omega_{b}$     $\omega_{dm}$     $\omega_{de}$

when $z\in[0,2]$     $73.52 \pm 0.24$ km/s/Mpc     $0.271 \pm 0.012$     $0.0036 ^{+0.0038}_{-0.0025}$      $-0.0505 \pm 0.0097$     $ -1.049 \pm 0.044$

when $z\in[0,2.5]$     $73.47 \pm 0.23$ km/s/Mpc     $0.268 \pm 0.011$     $0.0036 ^{+0.0041}_{-0.0022}$     $-0.0354 \pm 0.0093$     $-1.091 \pm 0.039$

when $z\in[0,5]$      $73.65 \pm 0.15$ km/s/Mpc     $0.254 \pm 0.008$     $0.0069 \pm 0.0045$     $-0.0752 \pm 0.0047$     $ -1.064 \pm 0.009$

Best fit values and $1\sigma$ errors estimated for the X (Dark Energy with $\omega_{de} \neq -1$) + Baryonic Matter + Dark Matter (with $\omega_{dm} \neq 0$) model, when $z \in [0,2]$, $z \in [0,2.5]$ and $z \in [0,5]$, respectively. The results have been obtained from the Bayesian (Probabilistic) Machine Learning approach, where the generative based process is based on Eqs. (<ref>) and (<ref>), using $H_{0} \in [64.0, 80.0]$,  $\Omega_{dm} \in [0.23, 0.4]$, $\Omega_{b} \in [0.0, 0.1]$, $\omega_{dm} \in [-0.1, 0.1]$, and $\omega_{de} \in [-1.5, -0.4]$ as flat priors. The analysis is based on 10 chains and in each chain, 10,000 "observational" data-sets from the model have been simulated or generated.

Our second model is the X (Dark Energy with $\omega_{de} \neq -1$) + Baryonic Matter + Dark Matter (with $\omega_{dm} \neq 0$) model where $\omega_{de}$ needs not be fixed to be the cosmological constant. Instead, it is a barotropic dark fluid with the constant equation of state parameter, which should be fitted. A detailed analysis of the model by using a Bayesian Learning approach shows that:

* The best fit values for the model's free parameters with $1\sigma$ error is $H_{0} = 73.52 \pm 0.24$ km/s/Mpc, $\Omega_{dm} = 0.271 \pm 0.012$, $\Omega_{b} = 0.0036 ^{+0.0038}_{-0.0025}$, $\omega_{dm} = -0.0505 \pm 0.0097$ and $\omega_{de} = -1.049 \pm 0.044$, when $z \in [0,2]$. The contour map is depicted in Fig. (<ref>), in purple colour.

* As a second case, when $z \in [0,2.5]$, the most likely best fit values with $1\sigma$ error are: $H_{0} = 73.47 \pm 0.23$ km/s/Mpc, $ \Omega_{dm} = 0.268 \pm 0.011$, $\Omega_{b} = 0.0036 ^{+0.0041}_{-0.0022}$, $\omega_{dm} = -0.0354 \pm 0.0093$, and $\omega_{de} = -1.091 \pm 0.039$. The contour map is given in Fig. (<ref>), in orange colour.

* Eventually, when $z \in [0,5]$, the most likely best fit values with $1\sigma$ error are: $H_{0} = 73.65 \pm 0.15$ km/s/Mpc, $ \Omega_{dm} = 0.254 \pm 0.008$, $\Omega_{b} = 0.0069 \pm 0.0045$, $\omega_{dm} = -0.0752 \pm 0.0047$, and $\omega_{de} = -1.064 \pm 0.009$. The contour map is given in Fig. (<ref>), in red colour.

\begin{array}{cccc}
\includegraphics[width=120 mm]{model_XCDM.pdf}
\end{array}$

Contour maps of the model given by Eqs. (<ref>) and (<ref>), for $z \in [0,2]$, $z \in [0,2.5]$, and $z \in [0,5]$, respectively. The best fit values of the model's parameters have been found to be $H_{0} = 73.52 \pm 0.24$ km/s/Mpc, $\Omega_{dm} = 0.271 \pm 0.012$, $\Omega_{b} = 0.0036 ^{+0.0038}_{-0.0025}$, $\omega_{dm} = -0.0505 \pm 0.0097$, and $\omega_{de} = -1.049 \pm 0.044$, when $z \in [0,2]$. While, when $z \in [0,2.5]$, the best fit values of the parameters have been found to be  $H_{0} = 73.47 \pm 0.23$ km/s/Mpc, $ \Omega_{dm} = 0.268 \pm 0.011$, $\Omega_{b} = 0.0036 ^{+0.0041}_{-0.0022}$, $\omega_{dm} = -0.0354 \pm 0.0093$, and $\omega_{de} = -1.091 \pm 0.039$. Finally, the best fit values of the model parameters have been found to be $H_{0} = 73.65 \pm 0.15$ km/s/Mpc, $ \Omega_{dm} = 0.254 \pm 0.008$, $\Omega_{b} = 0.0069 \pm 0.0045$, $\omega_{dm} = -0.0752 \pm 0.0047$ and $\omega_{de} = -1.064 \pm 0.009$, when $z \in [0,5]$. In all three cases, $H_{0} \in [64.0, 80.0]$,  $\Omega_{dm} \in [0.23, 0.4]$, $\Omega_{b} \in [0.0, 0.1]$, $\omega_{dm} \in [-0.1, 0.1]$ and $\omega_{de} \in [-1.5, -0.4]$ flat priors have been imposed during the generative process used to generate the "observational" data. The analysis is based on 10 chains and in each chain, 10,000 "observational" data-sets from the model have been simulated or generated.

We conclude from the results above that the $H_{0}$ tension can be solved and that there is a clear deviation from cold dark matter. Moreover, we see that this deviation is stronger than in previous case. Even more, the solution does not require any non-gravitational interaction between dark energy and dark matter. On the other hand, we observe that the deviation from cold dark matter has a clear redshift dependence and, moreover, if the expansion rate data for $z \in [0,5]$ would become available, then we should expect a very strong deviation. The results of the forecast for $z \in[0,5]$ gives us the following $H_{0} = 73.65 \pm 0.15$ km/s/Mpc, $ \Omega_{dm} = 0.254 \pm 0.008$, $\Omega_{b} = 0.0069 \pm 0.0045$, $\omega_{dm} = -0.0752 \pm 0.0047$ and $\omega_{de} = -1.064 \pm 0.009$ constraints. We notice that, for the higher redshift expansion rate data, the mean of $\Omega_{dm}$ will decrease significantly. For instance, if it is $\Omega_{dm} = 0.271$ for $z \in [0,2]$, then it will be $\Omega_{dm} = 0.268$ for $z \in [0,2.5]$. While we have obtained  $\Omega_{dm} = 0.254$, when $z\in [0,5]$.  On the other hand, the constraints on $\omega_{de}$ show that an increase in $\omega_{dm}$ will mean a decrease in $\omega_{de}$. Additionally, the  constraints obtained on $\omega_{de}$ in all three redshift ranges are in good agreement with existing observational data reported in the literature. Finally, we observe that the expansion rate data for $z\in [0,5]$ will indicate a significant change in $\Omega_{b}$ (see Table <ref>).

To summarize our results, let us briefly mention where we started and what we have achieved. We have an $H_{0}$ tension and see that interacting dark energy models can solve it. If we have a closer look at introducing non-gravitational interaction between dark energy and dark matter, then it becomes clear that it affects the equation of state parameters of two energy sources. On the other hand, we have one extra degree of freedom, making the model softer for observational data. The question is why this should work and be useful for the $H_{0}$ tension problem? Can we have another look at the problem from a different perspective and figure out how it works? Here we worry only about the non-gravitational interaction since it is not clear why it should be used and how it appears in our Universe. Of course, we could argue that the observational data supports this idea. However, we should remember that it is true and workable when we have a specific dark energy model and cold dark matter, which is actually the only way to gauge it and does not give answers.

Now, what we did in this paper. In particular, taking into account that the non-gravitational interaction affects the equation of state parameters of dark energy and dark matter models, we have applied Bayesian (Probabilistic) Machine Learning to study two cosmological models where, in each case, $\omega_{dm}$ is the free parameter that should be fitted. Surprisingly we have learned that $\omega_{dm} \neq 0$, both when dark energy is the usual cosmological constant and when it is given by another energy source, with $\omega_{de} \neq -1$. In both cases, we have learned the constraints on the free parameters and concluded that the $H_{0}$ tension problem is strongly dependent on the non-coldness of dark matter. Therefore, making dark matter colder for a given dark energy model without using non-gravitational interaction allows to solve (or at least alleviate) the $H_{0}$ tension problem.

Eventually, we would like to mention that the obtained constraints indicate how the interactions should be built. In particular, starting from the models not containing non-gravitational interaction, we found that, on the different cosmological scales, there is a clear and increasing deviation from the cold dark matter paradigm, which in the light of the $H_{0}$ tension problem could be confused to be an interaction between dark energy and dark matter. The graphic behavior of the expansion rate $H(z)$ and the deceleration parameter $q(z)$, corresponding to this case, can be found in Fig. (<ref>). Similarly to the first case, here we see that very good learning has been achieved. Moreover, we see that the learned deviation for $\omega_{dm}$ is not only able to solve the $H_{0}$ tension, but it can also explain the BOSS experiment result of the expansion rate data at $z=2.34$ (the purple and orange curves in the left-hand side plot of Fig. (<ref>)). However, the forecast results for $z\in [0,5]$ hint that the high redshift expansion rate data can be in tension with the learned $\omega_{dm}$ deviation and could be not sufficient to explain the BOSS experiment result. Here, we should mention that this learned property can be used to impose the lowest bound on the $\omega_{dm}$ deviation. This is easy to see from a closer examination of the learned constraints presented in Table <ref>. Moreover, from Table <ref>, we see that we have learned very tight constraints for the model's free parameters, except $\Omega_{b}$.

\begin{array}{cccc}
\includegraphics[width=80 mm]{H_Model_2.jpeg}&&
\includegraphics[width=80 mm]{q_Model_2.jpeg}\\
\end{array}$

The graphic behavior of the Hubble parameter, in comparison with known $H(z)$ data, is presented on the left-hand side. The purple curve is a plot of the Hubble parameter for the best fit values of the model's parameters, when $z \in [0,2]$; the dashed orange curve corresponds to the case when $z \in [0,2.5]$; while the dashed red curve is for the case when $z \in [0,5]$. The red dots correspond to the known observational $H(z)$ data, and are the same as in Table <ref>. On the right-hand side, the graphic behavior of the deceleration parameter $q(z)$ is plotted. The model is given by Eqs. (<ref>) and  (<ref>). In both cases, only the best fit values for the model's free parameters obtained by the Bayesian Learning approach and presented in Table <ref> have been used.

§ CONCLUSIONS

Bayesian (Probabilistic) Machine Learning is one of the approaches used in modern science. It allows the use of the model itself to learn about its specific properties. It is a generative-based approach allowing us to explore the regimes in which we do not have physical data or where collecting them is a very costly process. It has been shown already that this procedure can be beneficial for other learning architectures by significantly reducing the associated costs. In this regard, the method can be very useful for predictive purposes too. In this paper, we use Bayesian (Probabilistic) Machine Learning to learn if and why interacting dark energy models should be useful in view of the $H_{0}$ tension problem. A closer look at interacting dark energy models indicates that the interaction does affect the equation of state parameters describing the two energy sources. However, why this interaction between dark energy and dark matter should exist? Could it be possible that we did not use proper parametrizations to describe the dark energy and dark matter components and that the interaction itself does only serve to correct this? Is this the reason why it seems to be quite useful to solve different cosmological problems?

To better understand the situation, we decided to learn how precise our assumptions constrain the equation of state parameters of dark energy and cold dark matter. We started assuming that the dark matter $\omega_{dm}$ equation of state parameter is one of the free parameters to be fitted in the context of the $\Lambda$ + Baryonic Matter + Dark Matter (with $\omega_{dm} \neq 0$) and  X (Dark Energy with $\omega_{de} \neq -1$) + Baryonic Matter + Dark Matter (with $\omega_{dm} \neq 0$) models. The constraints have been obtained for three growing redshift ranges, $z\in [0,2]$, $z\in [0,2.5]$ and $z\in [0,5]$, respectively. The first two redshift ranges are chosen to validate our fit results since we have very good quality observations of the expansion rate $H(z)$ in these cases. On the other hand, considering the extended redshift range, $z\in [0,5]$, is for forecasting reasons. The full validation of the results will have to wait for available observations of higher-redshift $H(z)$ data.

Surprisingly, we learned that $\omega_{dm} \neq 0$, both when dark energy is the cosmological constant and an energy source with $\omega_{de} \neq -1$. We clearly see a deviation from the $\Lambda$CDM standard model. The learning results based on the generated expansion rate data from the background dynamics of each model show that, on the cosmological scales, there is a deviation from the cold dark matter hypothesis. The results obtained explicitly demonstrate that the deviation from cold dark matter can effectively be a solution for the $H_{0}$ tension problem. In particular, we learned that, by only making dark matter colder for a given dark energy model different from the cosmological constant without using a non-gravitational interaction, one might get a solution to the $H_{0}$ tension problem. Moreover, our results hint at why the parameters of the interaction term should be small, as often discussed in the literature.

On the other hand, we saw that the learned deviation for $\omega_{dm}$ for the $\Lambda$ + Baryonic Matter + Dark Matter (with $\omega_{dm} \neq 0$) can not solve the $H_{0}$ tension, but it is able to explain the BOSS experiment result of the expansion rate data at $z=2.34$. It should be mentioned that the forecast results for $z\in [0,5]$ allow us to conclude that the high redshift expansion rate data can validate that the learned $\omega_{dm}$ deviation still may be useful to explain the BOSS experiment result. Contrary to the first case, the analysis of the X (Dark Energy with $\omega_{de} \neq -1$) + Baryonic Matter + Dark Matter (with $\omega_{dm} \neq 0$) model shows that the learned deviation for $\omega_{dm}$ does the job. However, the forecast results for $z\in [0,5]$ hint that the high-redshift expansion rate data could be in tension with the learned $\omega_{dm}$ deviation and that they might not be sufficient to explain the BOSS experiment result. Here, we should mention that this learned property can impose the lowest bound on the $\omega_{dm}$ deviation. This is easy to see, particularly from a closer examination of the learned constraints presented in Table <ref>. The results of Table <ref> and those of  Table <ref> indicate that we can learn very tight constraints on the free parameters, except for $\Omega_{b}$.

Finally, we should duly stress that, while playing with the equation of state parameter of dark matter, we have achieved some remarkable results indicating a clear deviation from the cold dark matter paradigm. The last result could be applied to solve the $H_{0}$ tension indicating new physics, which should be explored further. They could be applied to explain the BOSS experimental result, too. In the present work, we have restricted ourselves to two specific dark energy models only. However, other viable dark energy models are available,  and it will be most interesting to investigate them to learn more about the $\omega_{dm}$ deviation. Therefore, we plan to extend our work to consider more complicated cases involving other dark energy models. We also plan to study the structure formation problem, which is key to understanding how a deviation from cold dark matter will affect structure formation and stability. Of course, understanding to which extend the learned deviation can be useful to overcome various cosmological problems (as a replacement for involving interacting dark energy models) is a crucial issue that we plan to address in new situations. Progress on the mentioned issues will be reported in forthcoming papers.  To end, we would like to recall that the whole analysis presented in this paper has been done employing a PyMC3 python-based framework, where the expansion rate data has been generated and used during the learning process. Only after that, the learned constraints on the model's free parameters have been used to plot the graphical behaviour of the $H(z)$ parameter for each model to compare it with the available observational data.

§ ACKNOWLEDGEMENTS

This work has been partially supported by MCI/AEI/FEDER, European Union, projects FIS2016-76363-P and PID2019-104397GB-I00, by AGAUR (Catalan Government), project 2017-SGR-247, by the Polish National Science Centre (NCN) under the Grant Agreement 2020/37/B/ST2/02371 and the research activities co-financed by the funds granted under the Research Excellence Initiative of the University of Silesia in Katowice.