Fast Graph Fourier Transforms Based on Graph Symmetry and Bipartition

TL;DR

This paper introduces fast algorithms for graph Fourier transforms leveraging graph symmetry and bipartition, significantly reducing computational costs for symmetric and structured graphs in applications like image processing and spectral clustering.

Contribution

It develops novel methods using Haar units and graph topological symmetry to accelerate GFT computation, especially for symmetric and nearly regular graphs.

Findings

Significant reduction in computation costs demonstrated.

Effective for symmetric and structured graphs like line, cycle, and skeletal graphs.

Applicable in video compression and human action analysis.

Abstract

The graph Fourier transform (GFT) is an important tool for graph signal processing, with applications ranging from graph-based image processing to spectral clustering. However, unlike the discrete Fourier transform, the GFT typically does not have a fast algorithm. In this work, we develop new approaches to accelerate the GFT computation. In particular, we show that Haar units (Givens rotations with angle $\pi/4$) can be used to reduce GFT computation cost when the graph is bipartite or satisfies certain symmetry properties based on node pairing. We also propose a graph decomposition method based on graph topological symmetry, which allows us to identify and exploit butterfly structures in stages. This method is particularly useful for graphs that are nearly regular or have some specific structures, e.g., line graphs, cycle graphs, grid graphs, and human skeletal graphs. Though butterfly stages based on graph topological symmetry cannot be used for general graphs, they are useful in applications, including video compression and human action analysis, where symmetric graphs, such as symmetric line graphs and human skeletal graphs, are used. Our proposed fast GFT implementations are shown to reduce computation costs significantly, in terms of both number of operations and empirical runtimes.

Tables

Table 1. TABLE I: Definitions of symmetries for vectors, matrices, and graphs.

Subject	Terminology	Definition
Vector $𝐯$	Even symmetric	$𝐯=\mathrm{𝐉𝐯}$
	Odd symmetric	$𝐯=-\mathrm{𝐉𝐯}$
Matrix $𝐌$	Symmetric	$𝐌={𝐌}^{\top }$
	Centrosymmetric [35]	$𝐌={\mathrm{𝐉𝐌}}^{\top }𝐉$
	Bisymmetric [35]	$𝐌={𝐌}^{\top }={\mathrm{𝐉𝐌}}^{\top }𝐉$
Graph	$\varphi$-symmetric (Definition 2)	$w_{i,j}=w_{\varphi \left(i\right),\varphi \left(j\right)},\forall i,j$
$𝒢(𝒱,ℰ,𝐖)$

Table 2. TABLE II: Types of symmetric N × N 𝑁 𝑁 N\times N grids and their corresponding involutions. Indices k 𝑘 k and l 𝑙 l represent vertical and horizontal coordinates of grid nodes, as in Fig.  11 .

Symmetry type	Involution
Centrosymmetry	$\varphi \left((k,l)\right)=(N+1-l,N+1-k)$
UD-symmetry	$\varphi \left((k,l)\right)=(N+1-k,l)$
LR-symmetry	$\varphi \left((k,l)\right)=(k,N+1-l)$
Diagonal symmetry	$\varphi \left((k,l)\right)=(l,k)$
Anti-diagonal symmetry	$\varphi \left((k,l)\right)=(N+1-l,N+1-k)$

Table 3. TABLE III: Speed performance of proposed fast GFTs. The baseline for the runtime reduction rates is the matrix GFT implementation.

Topology	$n$	Number of Operations		Runtime
		Matrix ($\pm$/$\ast$)	Fast ($\pm$/$\ast$)	Reduction
Cycle	12	132$/$144	44$/$30	52.7%
	80	6320$/$6400	1224$/$1078	79.7%
6-conn. grid	16	240$/$256	80$/$80	53.7%
	64	4032$/$4096	1104$/$1072	68.5%
Z-shaped grid	16	240$/$256	128$/$112	41.5%
	64	4032$/$4096	2048$/$2048	45.0%
Skeleton	15	210$/$225	96$/$102	45.5%
	25	600$/$625	272$/$282	47.5%

Equations

$$\left\{\begin{array}[]{ll}\Theta_{pp}=\Theta_{qq}=\cos\theta,\\ \Theta_{qp}=-\Theta_{pq}=\sin\theta,\\ \Theta_{ii}=1,&i\neq p,\;q,\\ \Theta_{ij}=0,&\text{otherwise}.\end{array}\right.$$

$${\mathbf{J}}_n=\left(\begin{array}{cccc}& & & 1\\ & & 1& \\ & \text{udots}& & \\ 1& & & \end{array}\right).$$

$$\text{diag}(a_1,a_2,\dots ,a_k)=\left(\begin{array}{cccc}{a}_1& & & \\ & a_2& & \\ & & \ddots & \\ & & & a_k\end{array}\right),$$

$$\text{diag}({\mathbf{M}}_1,{\mathbf{M}}_2,\dots ,{\mathbf{M}}_k)=\left(\begin{array}{cccc}{\mathbf{M}}_1& & & \\ & {\mathbf{M}}_2& & \\ & & \ddots & \\ & & & {\mathbf{M}}_k\end{array}\right).$$

$$\mathbf{L}=\mathbf{D}-\mathbf{W}+\mathbf{S},$$

$$s_i=\sum _{j=1}^n{l}_{i,j},w_{i,j}=-l_{i,j}\text{ for }i\ne j,$$

$$l_{i,i}=s_i+\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^n{w}_{i,j},l_{i,j}=-w_{i,j}\text{ for }i\ne j.$$

$${\mathbf{f}}^{\top }\mathbf{L}\mathbf{f}=\sum _{(i,j)\in \mathcal{E}}{w}_{i,j}(f_i-f_j{)}^2+\sum _{k=1}^n{s}_k{f}_k^2$$

$${\mathbf{u}}_1=\underset{\Vert \mathbf{f}\Vert =1}{\mathrm{argmin}}{\mathbf{f}}^{\top }\mathbf{L}\mathbf{f},{\mathbf{u}}_k=\underset{\mathbf{f}\perp {\mathbf{u}}_1,\dots ,{\mathbf{u}}_{k-1},\Vert \mathbf{f}\Vert =1}{\mathrm{argmin}}{\mathbf{f}}^{\top }\mathbf{L}\mathbf{f}.$$

$$\left(\begin{array}{c}{\hat{x}}_1\\ {\hat{x}}_2\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)\cdot {\Theta }^{\top }(1,2,\pi /4)\cdot \left(\begin{array}{c}{x}_1\\ x_2\end{array}\right).$$

$${\mathbf{U}}_{{\mathcal{C}}_4}=\frac{1}{2}\left(\begin{array}{cccc}1& 1& 1& 1\\ 1& -1& 1& -1\\ 1& -1& -1& 1\\ 1& 1& -1& -1\end{array}\right).$$

$${\mathbf{B}}_{n,p}=\frac{1}{\sqrt{2}}\left(\begin{array}{ccc}{\mathbf{I}}_p& 0& {\mathbf{J}}_p\\ 0& \sqrt{2}{\mathbf{I}}_{n-2p}& 0\\ {\mathbf{J}}_p& 0& -{\mathbf{I}}_p\end{array}\right).$$

$$({\bf B}_{n,p}\cdot{\bf x})_{i}=\left\{\begin{array}[]{ll}\frac{1}{\sqrt{2}}(x_{i}+x_{n+1-i}),&i=1,\dots,p\\ x_{i},&i=p+1,\dots,n-p\\ \frac{1}{\sqrt{2}}(-x_{i}+x_{n+1-i}),&i=n-p+1,\dots,n\end{array}\right.$$

$${\mathbf{U}}^{\top }\mathbf{x}={\mathbf{B}}_{n,p}\cdot diag({\mathbf{E}}^{\top },{\mathbf{F}}^{\top })\cdot \mathbf{x},$$

$$\mathbf{U}$$

$$=\frac{1}{\sqrt{2}}\left(\begin{array}{ccc}{\mathbf{E}}_{11}& \sqrt{2}{\mathbf{E}}_{12}& {\mathbf{E}}_{11}{\mathbf{J}}_p\\ {\mathbf{E}}_{21}& \sqrt{2}{\mathbf{E}}_{22}& {\mathbf{E}}_{21}{\mathbf{J}}_p\\ \mathbf{F}{\mathbf{J}}_p& 0& -\mathbf{F}\end{array}\right),$$

$$diag(\mathbf{E},\mathbf{F})=\mathbf{U}\cdot {\mathbf{B}}_{n,p}.$$

$$\mathbf{U}={\mathbf{B}}_{n,n/2}\left(\begin{array}{cc}{\mathbf{U}}^{+}& 0\\ 0& {\mathbf{U}}^{-}\end{array}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}{\mathbf{U}}^{+}& \mathbf{J}{\mathbf{U}}^{-}\\ \mathbf{J}{\mathbf{U}}^{+}& -{\mathbf{U}}^{-}\end{array}\right),$$

$$\mathbf{G}$$

$$\mathbf{L}=\left(\begin{array}{ccc}{\mathbf{L}}_{XX}& {\mathbf{L}}_{XZ}& {\mathbf{L}}_{XY}\\ {\mathbf{L}}_{ZX}& {\mathbf{L}}_{ZZ}& {\mathbf{L}}_{ZY}\\ {\mathbf{L}}_{YX}& {\mathbf{L}}_{YZ}& {\mathbf{L}}_{YY}\end{array}\right),\mathbf{G}=\left(\begin{array}{ccc}{\mathbf{G}}_{XX}& {\mathbf{G}}_{XZ}& {\mathbf{G}}_{XY}\\ {\mathbf{G}}_{ZX}& {\mathbf{G}}_{ZZ}& {\mathbf{G}}_{ZY}\\ {\mathbf{G}}_{YX}& {\mathbf{G}}_{YZ}& {\mathbf{G}}_{YY}\end{array}\right).$$

$$\mathbf{U}={\mathbf{B}}_{n,p}\cdot diag({\mathbf{U}}^{+},{\mathbf{U}}^{-}),{\mathbf{U}}^{+}\in {\mathbb{O}}^{n-p},{\mathbf{U}}^{-}\in {\mathbb{O}}^p.$$

$${\mathbf{L}}_{YY}=\mathbf{J}{\mathbf{L}}_{XX}\mathbf{J},{\mathbf{L}}_{YX}=\mathbf{J}{\mathbf{L}}_{XY}\mathbf{J},{\mathbf{L}}_{ZY}={\mathbf{L}}_{ZX}\mathbf{J}.$$

$$\mathbf{U}=\frac{1}{\sqrt{2}}\left(\begin{array}{ccc}{\mathbf{I}}_p& 0& {\mathbf{J}}_p\\ 0& \sqrt{2}{\mathbf{I}}_{n-2p}& 0\\ {\mathbf{J}}_p& 0& -{\mathbf{I}}_p\end{array}\right)\left(\begin{array}{ccc}{\mathbf{U}}_{XX}^{+}& {\mathbf{U}}_{XZ}^{+}& 0\\ {\mathbf{U}}_{ZX}^{+}& {\mathbf{U}}_{ZZ}^{+}& 0\\ 0& 0& {\mathbf{U}}^{-}\end{array}\right)$$

$${\mathbf{G}}_{XY}$$

$${\mathbf{G}}_{ZY}$$

$$\mathbf{G}=diag({\mathbf{V}}_1,{\mathbf{V}}_2)\cdot diag({\lambda }_1,\dots ,{\lambda }_n)\cdot diag({\mathbf{V}}_1,{\mathbf{V}}_2{)}^{\top },$$

$$\mathbf{G}=\left(\begin{array}{ccc}{\mathbf{L}}_{XX}+{\mathbf{L}}_{XY}\mathbf{J}& \sqrt{2}{\mathbf{L}}_{XZ}& 0\\ \sqrt{2}{\mathbf{L}}_{XZ}^{\top }& {\mathbf{L}}_{ZZ}& 0\\ 0& 0& {\mathbf{L}}_{YY}-\mathbf{J}{\mathbf{L}}_{XY}\end{array}\right),$$

$${\mathbf{L}}^{+}:=\left(\begin{array}{cc}{\mathbf{L}}_{XX}+{\mathbf{L}}_{XY}\mathbf{J}& \sqrt{2}{\mathbf{L}}_{XZ}\\ \sqrt{2}{\mathbf{L}}_{XZ}^{\top }& {\mathbf{L}}_{ZZ}\end{array}\right),{\mathbf{L}}^{-}:={\mathbf{L}}_{YY}-\mathbf{J}{\mathbf{L}}_{XY}.$$

$$w_{i,j}=w_{n+1-i,n+1-j},\forall i\in \mathcal{V},j\in \mathcal{V}.$$

$$p_{\varphi }:=\frac{1}{2}\times \left|\{i\in \mathcal{V}:i\ne \varphi (i)\}\right|.$$

$$({\bf B}_{\phi})_{i,j}=\left\{\begin{array}[]{ll}1/\sqrt{2},&i=j\in{\cal V}_{X}\\ -1/\sqrt{2},&i=j\in{\cal V}_{Y}\\ 1,&i=j\in{\cal V}_{Z}\\ 1/\sqrt{2},&i\in{\cal V}_{X},\;j=\phi(i)\in{\cal V}_{Y}\\ 1/\sqrt{2},&i\in{\cal V}_{Y},\;j=\phi(i)\in{\cal V}_{X}\\ 0,&\text{otherwise}\end{array}\right.$$

$$\displaystyle w_{i,j}^{+}=\left\{\begin{array}[]{ll}w_{i,j}+w_{i,\phi(j)},&\text{ if }i\in{\cal V}_{X},\;j\in{\cal V}_{X}\\ \sqrt{2}w_{i,j},&\text{ if }i\in{\cal V}_{X},\;j\in{\cal V}_{Z}\text{ or }i\in{\cal V}_{Z},\;j\in{\cal V}_{X}\\ w_{i,j},&\text{ if }i\in{\cal V}_{Z},\;j\in{\cal V}_{Z},\end{array}\right.$$

Full text

I Introduction

Graph signal processing (GSP) [1, 2, 3] is a framework that extends signal processing tools to data lying on irregular domains. In GSP, data points are represented as nodes in a graph, and relations between data points are captured by the graph edges. Data associated to the nodes is called a graph signal. Conventional signal processing tools such as the Fourier transform and filtering can be extended to signals defined on graphs, providing applications in sensor networks [4], image and video processing [5], and machine learning [6].

As an extension of the discrete Fourier transform (DFT) to graph signals, the graph Fourier transform (GFT) is a fundamental tool in GSP. There are several definitions of the GFT [2, 1, 7], depending on whether the graph is directed, which graph shift operator is used (e.g. adjacency matrix or Laplacian matrix), and how the graph signal energy is defined. Following the definition in [2], the GFT basis functions are defined as the eigenvectors of the graph Laplacian matrix, and their associated frequencies are the corresponding eigenvalues. The GFT coefficients of a given signal ${\bf x}$ can be obtained by projecting ${\bf x}$ onto the GFT basis functions. Those coefficients corresponding to smaller eigenvalues (lower frequency) reflect the energy of signal components with smaller variation on the graph. The GFT has a wide range of applications. First, based on the GFT and its frequency interpretation, graph spectral filters [8] can be defined by multiplying the GFT coefficients by the frequency response in the graph spectral domain, leading to applications such as image denoising and edge-preserving smoothing [9, 10]. Second, when the graph signal is modeled by a Gaussian Markov random field (GMRF) [11], the corresponding GFT can be regarded as the optimal decorrelating transform for that class of signals. Based on this fact, the GFT has been applied to image and video compression [12, 13, 14]. Third, in machine learning, when relations between data points are modeled by a graph, the GFT can be used for data clustering [15] and dimensionality reduction for classification [16, 17].

Unlike the DFT, which can be implemented with the well-known fast Fourier transform (FFT) algorithm [18], in general there are no fast algorithms to compute GFTs. DFT basis functions are always even or odd symmetric, which can be exploited to obtain fast algorithms. In contrast, arbitrary graph topologies do not always lead to Laplacian eigenvectors with such symmetry properties. Lack of fast algorithms is a significant drawback for GSP approaches, particularly when the GFT needs to be applied repeatedly. This has led researchers to investigate techniques for fast GFT computation (see Section I-A for a review of recent work). In particular, Magoarou et. al. have proposed a series of approaches for fast GFTs [19, 20, 21], which use optimization techniques to approximate a GFT by a fast transform constructed with a series of of parallel Givens rotations.

Our work is motivated by noting that exact fast GFTs are available for certain graphs with particular structures. For example, the discrete cosine transform (DCT) is known to be the GFT of a line graph with uniform edge weights [22], and it has well-known fast algorithms [23, 24]. Another example is a butterfly structured implementation for Type-4 DST [25], whose corresponding graph is a line graph with uniform weights with an added self-loop in the first node. Because of the availability of fast algorithms, DCT and Type-4 DST have been adopted in codecs such as HEVC [26] and AV1 [27]. Motivated by these fast algorithms, in this paper our goal is to explore more general classes of graphs with fast GFTs. In our preliminary work [28, 29], we have introduced two classes of graphs, symmetric line and grid graphs, whose GFTs have a butterfly stage for fast implementation. This work extends and generalizes the results of [28, 29] to graphs that are bipartite or have more general symmetry properties. A more detailed outline of this paper will be presented in Section I-B.

I-A Related Work

An $n$ dimensional Givens rotation [30], commonly referred to as a butterfly [25, 20, 21], is a linear transformation that applies a rotation of angle $\theta$ to two coordinates, denoted as $p$ and $q$. Its associated matrix \hbox{\boldmath\Theta}(p,q,\theta) has the form:

[TABLE]

A system using layers of parallel Givens rotations (e.g., Fig. 1), can be used to design a fast approximate transform. In particular, each Givens rotation can be implemented using three lifting steps [31], which further reduces the number of operations involved.

Recently, several papers have focused on GFT-specific speedup techniques. The work in [19] uses a gradient-descent-based optimization approach to approximate the GFT matrix by a product of sparse matrices, while [20] refines this method such that the resulting transform matrix can approximately diagonalize the graph Laplacian. In [21], a truncated Jacobi algorithm was introduced for picking the Givens rotations used in the approximate fast GFT, leading to an implementation with the structure shown in Fig. 1. This approach was further analyzed in [32], which demonstrates that more Givens rotations are required to approximate the Laplacian eigenvectors whose corresponding eigenvalues are close. Although these methods [19, 20, 21] are able to find approximate fast GFTs, they do so without taking advantage of structural properties of the original graph.

I-B Contributions

The relation of topological properties of graphs, such as bipartition, repeated subgraphs, symmetry, and uniformity of weights, to the structure of the GFT bases is an important topic in GSP. In this work, we show that for graphs with certain symmetry or bipartition properties, exact and fast GFTs based on Haar units (butterflies with rotation angle $\pi/4$) can be designed. We propose divide-and-conquer fast GFT algorithms for symmetric graphs and demonstrate that the resulting fast GFTs lead to significant complexity reduction, potentially beneficial in hardware implementation or in scenarios where the graph is fixed and the corresponding GFT is applied multiple times. We show that graphs for which such Haar-unit-based fast GFTs can be developed are useful in applications such as video coding and human activity analysis.

Unlike fast approximate GFTs [19, 20, 21], our fast GFTs are based on graph topological properties, and are exact. Experimental results show that as long as the desired graph symmetry property is available, our fast GFTs can provide outperform the approach in [21] in terms of speed. With respect to our earlier work [28, 29], the main novelties of this paper are: 1) we define a notion of graph symmetry that gives rise to butterfly implementation of the GFT, and show that the results in [28, 29] are particular cases within this general framework; 2) we introduce, in addition to line and grid graphs, more examples of graphs with fast GFTs, such as star graphs, cycle graphs, and skeleton graphs; 3) we provide more comprehensive results, including experimental runtimes and comparisons with existing approaches.

The rest of this paper is organized as follows. Section II introduces notation and basic graph signal processing concepts. In Section III we derive the algebraic conditions for a GFT to have a left or right butterfly stage. In Section IV we define graph symmetry based on node pairings, and propose a graph decomposition method for designing fast GFTs. In Section V we show several examples of fast GFTs based on the proposed method, and highlight some applications of the derived fast GFTs. Section VI provides experimental results to demonstrate the runtime reduction provided by the fast GFTs. Finally, Section VII concludes this paper.

II Preliminaries

II-A Notations and Conventions

We use bold symbols to denote vectors and matrices. The $n\times n$ identity matrix is denoted by ${\bf I}_{n}$. The $n\times n$ order-reversal permutation matrix is denoted by

[TABLE]

When ${\bf J}_{n}$ right (resp. left) multiplies another matrix, it flips this matrix left to right (resp. up to down). In (2) and in what follows, the entries not included in the matrix are meant to be zero. The subscripts of ${\bf I}$ and ${\bf J}$ matrices indicate their sizes, and may be omitted for brevity. For scalars $a_{1},\dots,a_{k}$ and square matrices ${\bf M}_{i}$ with arbitrary sizes, we denote diagonal and block diagonal matrices in compact notation as

[TABLE]

Finally, the set of $n$-dimensional real-valued vectors is denoted as $\mathbb{R}^{n}$, and the set of $n\times n$ orthogonal matrices (with columns normalized to have unit norms) is denoted as $\mathbb{O}^{n}$.

II-B Graph Fourier Transform

In this paper, we focus on undirected graphs111We leave the extension to directed graphs for future work.. Let ${\bf x}$ be a length-$n$ graph signal associated to an undirected graph ${\cal G}({\cal V},{\cal E},{\bf W})$. In particular, there are $n$ nodes in the vertex set ${\cal V}$, each corresponding to an element of ${\bf x}$. Each edge $e_{ij}\in{\cal E}$ describes the inter-sample relation between nodes $i$ and $j$. ${\bf W}$ is the weighted adjacency matrix, whose $(i,j)$ entry, $w_{i,j}$, is the weight of the edge between nodes $i$ and $j$, and $s_{i}:=w_{ii}$ is the weight of the self-loop on node $i$. The (unnormalized) graph Laplacian matrix of ${\cal G}$ is:

[TABLE]

where ${\bf S}=\text{diag}(s_{1},\dots,s_{n})$ is the diagonal self-loop matrix, and the degree matrix ${\bf D}=\text{diag}(d_{1},\dots,d_{n})$ is a diagonal matrix with $d_{i}=\sum_{j=1}^{n}w_{i,j}$. A graph is bipartite if its vertices can be divided into two disjoint sets (or, two parts) ${\cal S}_{1}$ and ${\cal S}_{2}$ such that every edge connects a vertex in ${\cal S}_{1}$ and one in ${\cal S}_{2}$.

By definition of graph Laplacian (3), one can express the self-loop and edge weights in terms of entries of the Laplacian matrix ${\bf L}=(l_{i,j})_{i,j}$, and vice versa:

[TABLE]

The Graph Fourier Transform (GFT), also known as Graph-Based Transform (GBT), is obtained from the eigen-decomposition of the graph Laplacian matrix, {\bf L}={\bf U}\hbox{\boldmath\Lambda}{\bf U}^{\top}, where ${\bf U}$ is the matrix of eigenvectors and $\Lambda$ is the diagonal matrix of eigenvalues. The $i$-th coefficient of GFT of a graph signal ${\bf x}$ is defined as the projection of ${\bf x}$ onto ${\bf u}_{i}$, the $i$-th column of ${\bf U}$. In some applications, such as spectral clustering [15], it may be beneficial to use a GFT defined on the eigenvectors of the symmetric normalized Laplacian ${\bm{\mathcal{L}}}={\bf D}^{-1/2}{\bf L}{\bf D}^{-1/2}$. In what follows, we use GFTs associated to the unnormalized Laplacian matrix, unless stated otherwise.

GFT coefficients provide a frequency representation of the given signal, since GFT basis functions associated to lower (resp. higher) eigenvalues represent lower (resp. higher) variation on the graph. To see this, we note that the Laplacian quadratic form

[TABLE]

measures the variation of signal ${\bf f}$ on the graph. Since $w_{i,j}$ and $s_{k}$ are non-negative, ${\bf L}$ is positive semi-definite and thus ${\bf f}^{\top}{\bf L}{\bf f}$ is always non-negative. The eigenvectors of ${\bf L}$ are the solutions to

[TABLE]

Thus, eigenvectors ${\bf u}_{1}$, $\dots$, ${\bf u}_{n}$ form an orthogonal basis with functions having lower to higher variations on the graph. The quantities of their corresponding variations are given by the associated eigenvalues $\lambda_{1}$, $\dots$, $\lambda_{n}$, which are also called graph frequencies.

III Algebraic Conditions for Haar Units in GFTs

In this paper, we improve computation efficiency by using the elementary operation in Fig. 2, which is equivalent to a 2$\times$2 Haar transform. This operation is a Givens rotation with angle $\pi/4$, followed by a sign flip. In particular, in Fig. 2,

[TABLE]

We refer to the operator in Fig. 2 as Haar unit, as opposed to general Givens rotations, which are often referred to as “butterflies” [25, 33, 21]. We say that a butterfly stage is a stage in a transform diagram with several parallel Givens rotations or Haar units. For example, in Fig. 3(a), we call the stage that produces $y_{i}$ from $x_{i}$ a butterfly stage, and the operator that produces $y_{1}$ and $y_{8}$ from $x_{1}$ and $x_{8}$ a Haar unit of this butterfly stage. Note that the factor $1/\sqrt{2}$ of the Haar unit can usually be absorbed into other stages of the transform computation (see Fig. 3(a) as an example, where the factor $1/\sqrt{2}$ is merged into the later stage). Thus, a Haar unit typically requires an addition and a subtraction only.

We consider a divide and conquer framework based on stages of Haar units and parallel sub-transforms, as illustrated Fig. 3. For each Haar unit, we always assume that the two output variables, such as $y_{1}$ and $y_{8}$ in Fig. 3(a), will be inputs of different sub-transforms in the next stage. Otherwise, such a Haar unit, e.g., the one acting on $x_{1}$ and $x_{8}$, can be trivially absorbed into the next stage.

As a first example, we consider a 4-node cycle graph with no self-loops and unity edge weights as in Fig. 4(a). It has a GFT matrix

[TABLE]

Based on the structure of ${\bf U}_{{\cal C}_{4}}$, it can be seen that GFT can be implemented using two butterfly stages, as in Fig. 4(b). In what follows, we refer to first-stage Haar units (e.g., the one acting on $x_{1}$ and $x_{4}$) as left Haar units, and to those in the last stage (such as the one producing $\widehat{x_{1}}$ and $\widehat{x_{2}}$) as right Haar units. We will explore the conditions that allow a GFT to be factored into terms that include left and right Haar units, which will enable us to develop techniques for designing such fast GFTs. We will show that right Haar units are associated to bipartite graphs (Section III-A), while left Haar units are related to graph symmetries (Section III-B).

Note that the second and third columns of ${\bf U}_{{\cal C}_{4}}$ correspond to eigenvalue 2 (which has multiplicity 2). This means that the GFT basis is not unique, because we can obtain another orthogonal basis for the eigenspace corresponding to eigenvalue 2. An example of another basis for this GFT is the length-4 DFT, which has a well-known fast algorithm [18]. Despite this non-uniqueness, a fast algorithm for a particular GFT basis would still be useful: we can first apply it to obtain coefficients for this particular basis, then apply an $m$-dimensional rotation ($m\times m$ orthogonal transform) on those $m$ GFT coefficients associated to eigenvalues with a multiplicity $m>1$, to obtain coefficients associated to another GFT basis. For example, if we properly apply a rotation to $\widehat{x_{2}}$ and $\widehat{x_{3}}$ in Fig. 4(b), we can obtain the second and third DFT coefficients. In what follows, we study GFT implementations for which stages of Haar units are available. In cases where eigenvalues of high multiplicity are present, we favor the set of eigenvectors for the corresponding subspace that will lead to a more efficient implementation.

For a general graph with $n$ nodes, we define the following $n\times n$ orthogonal matrix to represent a stage of $p$ parallel Haar units (with $p\leq n/2$):

[TABLE]

Note that ${\bf B}_{n,p}^{\top}={\bf B}_{n,p}$, and when we multiply a vector ${\bf x}=(x_{1},\dots,x_{n})^{\top}$ by ${\bf B}_{n,p}$, we have

[TABLE]

For example, ${\bf B}_{8,4}$ are ${\bf B}_{8,3}$ are equivalent to the butterfly stages in Figs. 3(a) and (b), respectively, with a scaling constant $1/\sqrt{2}$. The factors $\sqrt{2}$ and $1/\sqrt{2}$ are included in (7) so that the columns of ${\bf B}_{n,p}$ have unit norms; in this way, when ${\bf U}$ is an orthogonal matrix and ${\bf U}={\bf B}_{n,p}\bar{{\bf U}}$, then $\bar{{\bf U}}$ is orthogonal as well, meaning that ${\bf U}$ can be factorized into a butterfly stage and another orthogonal transform.

III-A Conditions for Right Haar Units

Let an orthogonal transform ${\bf U}$ have a right butterfly stage with $p$ Haar units, and assume without loss of generality that the entries of input and output vectors are properly ordered. Then, in compact notation, the GFT of input ${\bf x}$ can be written as

[TABLE]

where ${\bf E}\in\mathbb{O}^{n-p}$ and ${\bf F}\in\mathbb{O}^{p}$. This means that

[TABLE]

where ${\bf E}_{11}$, ${\bf E}_{12}$, ${\bf E}_{21}$, and ${\bf E}_{22}$ are subblock components of ${\bf E}$. Recall that ${\bf J}$ flips a matrix left to right when right-multiplied. Thus, for $k=1,\dots,p$, if we denote the $k$-th column of ${\bf U}$ as ${\bf u}_{k}=({\bf e}_{k}^{\top},{\bf f}_{p-k+1}^{\top})^{\top}$, then the $(n-k+1)$-th column of ${\bf U}$ is $({\bf e}_{k}^{\top},-{\bf f}_{p-k+1}^{\top})^{\top}$. GFT matrices with this structure arise from k-regular bipartite graphs (k-RBGs):

Lemma 1** ([34]).**

Let ${\bf L}$ be the Laplacian of a k-RBG with ${\cal S}_{1}=\{1,\dots,n/2\}$ and ${\cal S}_{2}=\{n/2+1,\dots,n\}$. If ${\bf u}=({\bf u}_{1}^{\top},\;{\bf u}_{2}^{\top})^{\top}$ with ${\bf u}_{1},{\bf u}_{2}\in\mathbb{R}^{n/2}$ is an eigenvector of ${\bf L}$ with eigenvalue $\lambda$, then $\hat{{\bf u}}=({\bf u}_{1}^{\top},\;-{\bf u}_{2}^{\top})^{\top}$ is an eigenvector of ${\bf L}$ with eigenvalue $2k-\lambda$.

In this lemma, we fix $p=n-p=n/2$, which follows from the k-RBG topology. Although Lemma 1 was introduced in [34], where only unweighted graphs are considered, it can be trivially generalized for weighted graphs. Lemma 1 also leads to the following theorem.

Theorem 1**.**

Let ${\bf L}$ be the Laplacian of a k-RBG with ${\cal S}_{1}=\{1,\dots,n/2\}$ and ${\cal S}_{2}=\{n/2+1,\dots,n\}$, then there exists a GFT matrix that has the structure (III-A). Therefore, this GFT has a stage of right Haar units.

Proof: We prove this theorem by construction. Since the graph Laplacian is symmetric, it has $n$ linearly independent eigenvectors, which enables us to construct a set ${\cal H}$ of $n$ eigenvectors as follows.

• a.

${\cal H}\leftarrow\phi$.

• b.

Pick an eigenvector \hbox{\boldmath\mu}=(\hbox{\boldmath\mu}_{1}^{\top};\hbox{\boldmath\mu}_{2}^{\top})^{\top} of ${\bf L}$ that is not in the span of ${\cal H}$. Note that this vector is guaranteed to exist as long as ${\cal H}$ has less than $n$ elements.

• c.

If \hbox{\boldmath\mu}_{2}={\bf 0}, then let {\cal H}\leftarrow{\cal H}\cup\{\hbox{\boldmath\mu}\}. Otherwise, let \hbox{\boldmath\mu}^{\prime}=(\hbox{\boldmath\mu}_{1}^{\top};-\hbox{\boldmath\mu}_{2}^{\top})^{\top}, which, by Lemma 1, is also an eigenvector of ${\bf L}$ and does not belong to the span of ${\cal H}$. Then, we set {\cal H}\leftarrow{\cal H}\cup\{\hbox{\boldmath\mu},\hbox{\boldmath\mu}^{\prime}\}.

• d.

Repeat b. and c. until ${\cal H}$ has $n$ elements. ∎

Theorem 1 provides certain sufficient (but, in fact, not necessary) conditions for a GFT to have a stage of right Haar units, even though the GFT matrix may not be unique. The matrices ${\bf E}$ and ${\bf F}$ can be obtained by right multiplying the target GFT matrix ${\bf U}$ by ${\bf B}_{n,p}$:

[TABLE]

If we consider eigenvectors of a normalized Laplacian, a result similar to Theorem 1 can be derived.

Theorem 2**.**

Let ${\bm{\mathcal{L}}}$ be the normalized Laplacian of a bipartite graph with ${\cal S}_{1}=\{1,\dots,p\}$ and ${\cal S}_{2}=\{p+1,\dots,n\}$, then there exists a GFT matrix that has the structure (III-A).

We omit the proof for brevity. In this case, ${\cal G}$ only needs to be bipartite (rather than k-regular bipartite), and $p$ need not be $n/2$.

III-B Conditions for Left Haar Units

If $n$ is even and the GFT has a butterfly stage in the left with exactly $n/2$ Haar units as in Fig. 3(a), then

[TABLE]

where ${\bf U}^{+},{\bf U}^{-}\in\mathbb{O}^{n/2}$ denote non-zero block components that characterize two sub-transforms as in Fig. 3(a).222The symbols ${\bf U}^{+}$ and ${\bf U}^{-}$ are chosen for consistency with Haar units. They are used to denote sub-GFTs, as will become clear in Sec. IV-C. From the right hand side of (10) we see that each column ${\bf u}_{i}$ of ${\bf U}$ must be either even symmetric (i.e., ${\bf u}_{i}={\bf J}{\bf u}_{i}$) or odd symmetric (i.e., ${\bf u}_{i}=-{\bf J}{\bf u}_{i}$). In this case, the Laplacian must be centrosymmetric (symmetric around the center):

Lemma 2** ([35]).**

Let $n$ be even. An $n\times n$ matrix ${\bf Q}$ has a set of $n$ linearly independent eigenvectors that are even or odd symmetric if and only if ${\bf Q}$ is centrosymmetric, i.e., ${\bf Q}={\bf J}{\bf Q}^{\top}{\bf J}$.

Note that the Laplacian matrix ${\bf L}$ of an undirected graph is always symmetric (${\bf L}={\bf L}^{\top}$), but an additional centrosymmetry condition (${\bf L}={\bf J}{\bf L}^{\top}{\bf J}$) is required so that Lemma 2 holds. Such a matrix with both symmetries (${\bf L}={\bf L}^{\top}={\bf J}{\bf L}^{\top}{\bf J}$) is called bisymmetric, and its entries are symmetric around both diagonals. The various types of symmetries considered in this paper listed in Table I.

Lemma 2 states that for even $n$, a GFT can be factored to include $n/2$ left Haar units if and only if the associated Laplacian matrix is bisymmetric (with nodes properly ordered). We now generalize this result to the case when there are only $p<n/2$ Haar units in the first butterfly stage, and with a possibly odd $n$. Again, we assume without loss of generality that the graph nodes, input, and output variables are properly ordered (notations defined for general node ordering will be introduced in Section IV). We let ${\cal V}_{X}=\{1,\dots,p\}$, ${\cal V}_{Z}=\{p+1,\dots,n-p\}$, and ${\cal V}_{Y}=\{n-p+1,\dots,n\}$ be disjoint subsets of vertices. We define

[TABLE]

and denote the corresponding subblock components of ${\bf L}$ and ${\bf G}$ as

[TABLE]

Similar to (8) and (III-A), the GFT matrix with a first butterfly stage of $p$ Haar units has the form of

[TABLE]

Then the following lemma describes the conditions for ${\bf L}$ to have a GFT with $p$ left Haar units.

Lemma 3**.**

Let ${\bf L}$ be a graph Laplacian matrix, then there exists a GFT matrix ${\bf U}$ in the form of (13), i.e., the associated ${\bf G}_{YX}$, ${\bf G}_{YZ}$, ${\bf G}_{XY}$ and ${\bf G}_{ZY}$ are zero matrices, if and only if

[TABLE]

Note that when $Z$ is empty, i.e., $p=n/2$, (14) implies that ${\bf L}$ has to be centrosymmetric, as in Lemma 2.

Proof: If ${\bf U}$ is a GFT matrix satisfying (13), we denote the subblocks of ${\bf U}^{+}$ as ${\bf U}_{XX}^{+}$, ${\bf U}_{XZ}^{+}$, ${\bf U}_{ZX}^{+}$, and ${\bf U}_{ZZ}^{+}$, and rewrite (13) as

[TABLE]

Denote the matrix of eigenvalues of ${\bf L}$ as \hbox{\boldmath\Lambda}={\hbox{diag}}(\hbox{\boldmath\Lambda}_{X},\hbox{\boldmath\Lambda}_{Z},\hbox{\boldmath\Lambda}_{Y}) with subblock sizes $p$, $n-2p$, and $p$, respectively. Then, we can express each subblock of ${\bf L}$ by expanding {\bf L}={\bf U}\hbox{\boldmath\Lambda}{\bf U}^{\top} with (15), and we can trivially verify that (14) holds.

To show the converse, we assume that (14) holds. Expanding the right hand side of (11), we can express the subblocks of ${\bf G}$ in terms of those of ${\bf L}$. In particular,

[TABLE]

With these expressions, (14) implies that ${\bf G}_{XY}$, ${\bf G}_{ZY}$, and their transpose versions ${\bf G}_{YX}$, ${\bf G}_{YZ}$ are all zero. This means that ${\bf G}$ is block-diagonal, and thus has an eigendecomposition as

[TABLE]

where ${\bf V}_{1}\in\mathbb{O}^{n-p}$ and ${\bf V}_{2}\in\mathbb{O}^{p}$. It follows that ${\bf B}_{n,p}\cdot{\hbox{diag}}({\bf V}_{1},{\bf V}_{2})$ is an eigenmatrix of ${\bf L}={\bf B}_{n,p}\cdot{\bf G}\cdot{\bf B}_{n,p}^{\top}$ as in (13). ∎

Under the conditions of (14), ${\bf G}$ reduces to

[TABLE]

and ${\bf U}^{+}$ and ${\bf U}^{-}$ are respectively the eigenmatrices of

[TABLE]

A diagram with $n=8$, $p=3$ is shown in Fig. 3(b) as an example.

Note that the desired properties in ${\bf L}$ correspond to certain symmetry properties in the graph topology. If ${\cal V}_{Z}$ is empty, Lemma 3 implies that

[TABLE]

In this case, when we plot the nodes in order on a 1D line, we can identify an axis in the middle, around which all edges and self-loops are symmetric. An example of a graph whose Laplacian is bisymmetric is shown in Fig. 5(a).

More generally, if ${\cal V}_{Z}$ is nonempty, then the first two equations in (14) indicate that the sub-matrix of ${\bf L}$ associated to ${\cal V}_{X}$ and ${\cal V}_{Y}$ is bisymmetric. This means that ${\cal V}_{X}$ and ${\cal V}_{Y}$ contain vertices that are symmetric to each other. The third equation in (14) implies that when there is an edge connecting $k\in{\cal V}_{Z}$ and $i\in{\cal V}_{X}$, there must be an edge with the same weight connecting $k$ and $n+1-i\in{\cal V}_{Y}$ as well. An example of this type of graph is shown in Fig. 5(b), where ${\cal V}_{X}=\{1\}$, ${\cal V}_{Z}=\{2,3\}$, and ${\cal V}_{Y}=\{4\}$. Similar to Fig. 5(a), we can identify a symmetry around the middle, though nodes in ${\cal V}_{Z}$ are not paired with symmetric counterparts.

Based on the observations above, we see that left Haar units are available when the graph has symmetry properties related to Lemma 3. However, Lemma 3 assumes that the nodes had been ordered properly, so that the Laplacian has the required bisymmetric structure. In general, the node labels of a graph will not be such that this condition is automatically met, even if the graph is symmetric. For example, if a different node labeling is applied to the graph of Fig. 5(a), its corresponding Laplacian may not be bisymmetric anymore. In the next section we study methods to identify graph symmetries directly using node pairing functions, which will allow us to design fast GFT algorithms, regardless of how the nodes are initially labeled.

IV Fast GFTs Based on Graph Symmetry

In this section, we will characterize how Lemma 3 relates to the graph topology. In particular, we define the symmetry properties observed in Fig. 5 based on an involution (node-pairing function) in Section IV-A. Given an observed graph symmetry characterized by an involution, in Section IV-B we define the node sets ${\cal V}_{X}$, ${\cal V}_{Y}$, and ${\cal V}_{Z}$. In Section IV-C, we propose a graph decomposition approach for searching fast GFTs in stages.

IV-A Graph Symmetry Based on Node Pairing

The symmetries of Fig. 5 can be described in terms of complete (Fig. 5(a)) and incomplete (Fig. 5(b)) node pairings. Such pairings can be defined by bijective mappings that are their own inverses, namely, involutions:

Definition 1** ([36]).**

A permutation on a finite set ${\cal V}$ is called an involution if it is its own inverse, i.e., $\phi(\phi(i))=i$ for all $i\in{\cal V}$.

We will use them to identify graph symmetries.

Definition 2**.**

Let $\phi$ be an involution on the vertex set ${\cal V}$ of a graph ${\cal G}$, then ${\cal G}$ is $\phi$-symmetric if $w_{i,j}=w_{\phi(i),\phi(j)}$ for all $i\in{\cal V}$, $j\in{\cal V}$.

Note that in Definition 2, the required property has to hold also for $i=j$. That is, $s_{i}=w_{i,i}=w_{\phi(i),\phi(i)}=s_{\phi(i)}$, meaning that the self-loops on nodes $i$ and $\phi(i)$ are required to have the same weight. Also note that, among permutations, only involutions are valid for Definition 2, since the pairing functions that lead to the conditions in (14) can only be induced by involutions.333Note that graph symmetry can be defined differently in different contexts. In algebraic graph theory, graph symmetry is defined based on transitivity of vertices and edges [37]. In [38], a graph is called symmetric if there exists a non-identical permutation $\phi$ (not necessarily an involution) on the graph nodes that leaves the graph unaltered. These definitions are beyond the scope of this paper. When we refer to graph symmetry in this paper, we always assume an involution $\phi$ is specified such that Definition 2 holds.

Let ${\cal V}=\{1,\dots,n\}$ be the vertex set, and let us denote an involution $\phi$ as $\phi=(\phi(1),\phi(2),\dots,\phi(n))$. For example, the involutions corresponding to the symmetries of graphs in Figs. 5(a) and (b) are $\phi_{a}=(4,3,2,1)$ and $\phi_{b}=(4,2,3,1)$, respectively. We also denote the number of available Haar units for a given $\phi$ as

[TABLE]

IV-B Node Partitioning for Haar Units

Once we observe a graph symmetry and characterize it by an involution $\phi$, we can identify the nodes on the axis of symmetry, ${\cal V}_{Z}:=\{i\in{\cal V}:\;\phi(i)=i\}$, then partition the other nodes into two sets ${\cal V}_{X}$ and ${\cal V}_{Y}$ such that nodes in those sets belong to different sides of the symmetry axis. In this way, we can define an orthogonal matrix ${\bf B}_{\phi}$ as a permuted version of ${\bf B}_{n,p_{\phi}}$ based on $\phi$, ${\cal V}_{X}$, ${\cal V}_{Y}$, and ${\cal V}_{Z}$ in the following way:

[TABLE]

This means that ${\bf L}_{\phi}:={\bf B}_{\phi}^{\top}{\bf L}{\bf B}_{\phi}$ is a permuted version of (16), whose block diagonal structure gives the following theorem:

Theorem 3** (Block-diagonalization of Laplacian based on graph symmetry).**

Let the graph ${\cal G}$ with Laplacian ${\bf L}$ be $\phi$-symmetric. Then, $({\bf L}_{\phi})_{i,j}=({\bf L}_{\phi})_{j,i}=0$ if $i\in{\cal V}_{X}\cup{\cal V}_{Z}$ and $j\in{\cal V}_{Y}$.

While Theorem 3 is derived based on the unnormalized Laplacian ${\bf L}$, it holds for normalized Laplacian as well.

IV-C Main Approach–Decomposition of Symmetric Graphs

The block-diagonalization of (11) maps ${\bf L}$ to ${\bf G}$ via ${\bf B}_{n,p}$, with ${\bf G}$ in (16). Note that, from (4) and (5), we can draw a one-to-one correspondence between a matrix and a graph. In this way, ${\bf G}={\hbox{diag}}({\bf L}^{+},{\bf L}^{-})$ can be can regarded as the Laplacian of a graph with two connected components, denoted as ${\cal G}^{+}$ and ${\cal G}^{-}$, with Laplacians ${\bf L}^{+}$ and ${\bf L}^{-}$, vertex sets ${\cal V}^{+}:={\cal V}_{X}\cup{\cal V}_{Z}$ and ${\cal V}^{-}:={\cal V}_{Y}$, weight matrices ${\bf W}^{+}$ and ${\bf W}^{-}$ (possibly with negative weights), respectively. With this graph decomposition from ${\cal G}$ to ${\cal G}^{+}$ and ${\cal G}^{-}$, the GFT of ${\cal G}$ can be implemented by a butterfly stage ${\bf B}_{n,p}$, followed by the two sub-GFTs corresponding to ${\cal G}^{+}$ and ${\cal G}^{-}$. Explicitly considering the graphs resulting from this decomposition is useful because in some cases ${\bf L}^{+}$ and ${\bf L}^{-}$ may in turn have symmetry properties, which could be exploited to achieve additional reductions in complexity. Moreover, considering the transforms after the Haar units as GFTs could lead to better interpretations of the overall GFT.

Regarding ${\bf L}^{+}$ and ${\bf L}^{-}$ in (17) as graph Laplacians, we can use (4) and (5) to express self-loop and edge weights of ${\cal G}^{+}$ and ${\cal G}^{-}$ in terms of those of ${\cal G}$, as described in the following theorem.

Theorem 4**.**

If ${\cal G}$ is $\phi$-symmetric with node partitions ${\cal V}_{X}$, ${\cal V}_{Y}$, and ${\cal V}_{Z}$, then the weights of ${\cal G}^{+}$ (with vertex set ${\cal V}^{+}={\cal V}_{X}\cup{\cal V}_{Z}$) and ${\cal G}^{-}$ (with vertex set ${\cal V}^{-}={\cal V}_{Y}$) are given by

[TABLE]

Refer to Appendix -A for the proof. We use the toy examples of Figs. 6 and 7 (with ${\cal V}_{Z}=\emptyset$ and ${\cal V}_{Z}\neq\emptyset$, respectively) to illustrate the graph decomposition. Note that ${\cal G}^{+}$ and ${\cal G}^{-}$ may have negative weights even if ${\cal G}$ does not. For any signal ${\bf x}$, we denote the “sum” (low-pass) and “difference” (high-pass) outputs of Haar units as ${\bf x}^{+}$ and ${\bf x}^{-}$: $(({\bf x}^{+})^{\top},({\bf x}^{-})^{\top})^{\top}={\bf B}_{\phi}^{\top}{\bf x}$. For example, in Fig. 6, ${\bf x}^{+}=(y_{1},y_{2})^{\top}$ and ${\bf x}^{-}=(y_{3},y_{4})^{\top}$. In Fig. 7, ${\bf x}^{+}=(z_{1},z_{2},z_{3})^{\top}$ and ${\bf x}^{-}=(z_{4})$.

The graph construction of Theorem 4 creates two disconnected sub-graphs by removing all edges between ${\cal V}_{y}$ and ${\cal V}_{x}\cup{\cal V}_{z}$ and preserving all other edges, but changing some of the weights and adding self-loops. Three types of cases lead to one or two edges being removed:

Edges connecting two symmetric nodes $i\in{\cal V}_{X}$ and $\phi(i)\in{\cal V}_{Y}$. The edge with weight $b$ in Fig. 6(a) is an example of this case. These edges are removed and lead to self-loops with twice the original weight in ${\cal G}^{-}$ ($2b$ in this case). 2. 2.

Two symmetric edges: each connecting a node in ${\cal V}_{X}$ to a node in ${\cal V}_{Y}$. The two edges with weight $c$ in Fig. 6(a) are an example. These edges are removed, but lead to changes in two edge weights, with the weight of the edge in ${\cal G}^{+}$ increasing and that of the edge in ${\cal G}^{-}$ decreasing. Two self-loops are also added to the corresponding nodes in ${\cal G}^{-}$. 3. 3.

Two symmetric edges with a common node in ${\cal V}_{Z}$. Edges with weight $d$ in Fig. 7(a) belong to this case. This case results in a single edge being kept, with a modified edge weight and two self-loops in ${\cal G}^{+}$, and a self-loop in ${\cal G}^{-}$.

Note that, the signals ${\bf x}^{+}$ and ${\bf x}^{-}$ correspond to ${\cal G}^{+}$ and ${\cal G}^{-}$, and can be regarded as even and odd symmetric components of the original graph signal ${\bf x}$. An example associated to Fig. 7 is shown in Fig. 8. We can see that a graph signal can be decomposed into two components ${\bf x}_{\text{even}}$ and ${\bf x}_{\text{odd}}$, which correspond to ${\bf x}^{+}$ and ${\bf x}^{-}$, respectively, by

[TABLE]

In particular, ${\bf x}_{\text{even}}$ and ${\bf x}_{\text{odd}}$ have even and odd symmetries based on the node pairing, i.e. ${\bf x}_{\text{even}}(i)={\bf x}_{\text{even}}(\phi(i))$ and ${\bf x}_{\text{odd}}(i)=-{\bf x}_{\text{odd}}(\phi(i))$ for all $i\in{\cal V}$. This can be considered as a generalization of even and odd symmetric components decomposition for finite length time series. Components ${\bf x}_{\text{even}}$ and ${\bf x}_{\text{odd}}$ of the graph signal ${\bf x}$ can be regarded as intermediate results of the GFT coefficients.

The decomposition described in Theorem 4 enables us to search further stages of Haar units in the sub-GFTs ${\bf U}^{+}$ and ${\bf U}^{-}$ by inspecting their associated graphs ${\cal G}^{+}$ and ${\cal G}^{-}$. Once a symmetry based on an involution is found in ${\cal G}^{+}$ or ${\cal G}^{-}$, we can apply the decomposition again, and repeat until a symmetry property cannot be found anymore. Some examples will be provided in Section V.

V Examples and Applications

In practice, graphs with distinct weights on different edges or graphs learned from data without any topology constraints are not likely to have the desired bipartition and symmetry properties. However, bipartite and symmetric graph structures arise in graphs considered in certain fields. Examples of bipartite graphs include tree-structured graphs, whose GFTs are useful for designing wavelet transforms on graph [39]. Involution-based symmetries can be found in graphs with regular or partially regular topologies (e.g. line, cycle, and grid graphs), graphs that are symmetric by construction (e.g. human skeletal graphs), and uniformly weighted graphs. In what follows, we study several classes of graphs with these properties and discuss the search of involution in general graphs.

V-A Graphs with 2-Sparse Eigenvectors

The authors in [40] have studied the conditions for 2-sparse graph eigenvectors to exist:

Lemma 4** ([40]).**

A Laplacian has an eigenvector ${\bf u}$ with only two nonzero elements ${\bf u}(i)=1/\sqrt{2}$, ${\bf u}(j)=-1/\sqrt{2}$ if and only if

[TABLE]

In fact, the condition (22) is equivalent to having ${\cal G}$ $\phi$-symmetric, with $\phi(i)=j$, $\phi(j)=i$, and $\phi(k)=k$ for $k\neq i,j$. In this case, each of ${\cal V}_{X}=\{i\}$ and ${\cal V}_{Y}=\{j=\phi(i)\}$ has only one node, and ${\bf U}^{-}$ reduces to a one by one identity matrix. Examples of graphs satisfying Lemma 4 include uniformly weighted graphs with several types of topology: 1) star graph, 2) complete graph, and 3) a graph with a clique (a complete subgraph), where at least two nodes in the clique are not connected to any other nodes outside the clique.

V-B Symmetric Line Graphs

A Laplacian matrix associated to a line graph can be viewed as a precision matrix (inverse covariance matrix) of a first-order Gaussian Markov random field (GMRF), which can be used for modeling image and video pixels [41, 42]. One special case is the line graph with uniform weights, whose GFT is the well-known DCT.

If a line graph ${\cal G}_{l}$ is symmetric around the middle, then it is $\phi=(n,n-1,\dots,1)$-symmetric ${\cal G}_{l}$ and has a left butterfly stage. In our recent work [28], we consider inter-predicted residual blocks in video coding, where pixels in a residual block have nearly symmetric statistics around the middle. We model those blocks by a GMRF based on a symmetric line graph. The resulting GFT has a fast implementation and provides a coding gain as compared to the DCT.

V-C Steerable DFTs

The Laplacian of an $n$-node cycle graph ${\cal G}_{c}$ with unit weights is circulant and has non-unique GFTs since some of its eigenvalues have multiplicities greater than one. Due to the circulant structure, the DFT is one of the GFTs of ${\cal G}_{c}$. The family of all GFTs of ${\cal G}_{c}$ is called the steerable DFTs [43]. Note that, for any $n$, ${\cal G}$ is $\phi=(n,n-1,\dots,3,2,1)$-symmetric, which enables us to explore fast implementations for ${\cal G}_{c}$ other than the Cooley-Tukey fast Fourier transform (FFT) algorithm [18]. An example with $n=12$ is shown in Fig. 9, where two stages of Haar units are available, and some of the sub-GFTs after the first two stages can be further simplified. Unlike the conventional DFT, the derived GFT will have real operations only. In addition, while FFT cannot be easily applied when $n$ is a prime number, our method gives at least one left butterfly stage for any $n$. We also note that, for any steerable DFT with a length $n$ that is a multiple of 4, the GFTs of ${\cal G}_{c}^{++}$ and ${\cal G}_{c}^{-+}$ are Type-2 DCT and Type-4 DST, respectively. This means that those sub-GFTs can also be implemented using fast DCT and ADST algorithms [23, 24, 25]. To the best of our knowledge, fast implementations for steerable DFT other than the FFT algorithm have not been studied in the literature.

V-D Symmetric Grid Graphs

For the term “grid”, we refer to a graph with $n=N_{1}\times N_{2}$ nodes that correspond to integer positions in 2D Euclidean space. This means that each node can be associated to a 2D coordinate $(k,l)$ with $1\leq k\leq N_{1}$ and $1\leq l\leq N_{2}$. In practice, many grids that arise in applications such as image processing have highly regular topologies due to the structured data domain. For example, pixel data can be modeled by a 4-connected grid, where all nodes are connected to their 4 immediate neighbors only. When this grid is uniformly weighted, the 2D DCT is shown to be its GFT, and provides an optimal decorrelation of block data modeled by a 2D Gaussian Markov model [44]. In this paper, we focus on grids with nearly regular topologies (e.g., all internal nodes have the same number of neighbors) or particular symmetry properties.

GFTs on grids with arbitrary weights correspond to 2D non-separable transforms for pixel blocks, which can achieve a significant compression gain over the DCT [45]. In our recent work [29], we proposed speedup techniques for 2D grid-based transforms based on various types of grid symmetries: 1) centrosymmetry, 2) up down (UD-) symmetry, 3) left right (LR-) symmetry, 4) diagonal symmetry, 5) anti-diagonal, and 6) grids with multiple symmetry properties. These grid symmetries are defined based on different axes or point of symmetry, as shown in Fig. 11. While [29] applies different node re-ordering rules for different types of symmetric grid, in this work we can simply describe grid symmetries by the involutions shown in Table II, to provide a more straightforward derivation of fast GFTs on symmetric grids. For application, it has been shown in [46] that GFTs with exact or partial symmetry properties can enhance compression efficiency with respect to the DCT.

We provide two examples of fast GFTs on symmetric grids. The first example is shown in Fig. 10, with a 4$\times$4 grid ${\cal G}_{b}$ that is bi-diagonally symmetric (symmetric around both diagonals). We can first decompose ${\cal G}_{b}$ based on the diagonal symmetry into ${\cal G}_{b}^{+}$ and ${\cal G}_{b}^{-}$. Then, we observe that the symmetry around the anti-diagonal remains in ${\cal G}_{b}^{+}$ and ${\cal G}_{b}^{-}$, so further decomposition can be applied. As a result, the overall GFT has two butterfly stages of Haar units, and can be implemented using 4 sub-GFTs with length 6, 4, 4, and 2, as in Fig. 10. For the second example, we consider a grid graph in the coding framework proposed in [47]. This grid, which we refer to as z-shaped grid, is 4-connected grid with horizontal and anti-diagonal edges, as shown in the top-left of Fig. 12(a). We denote this grid as ${\cal G}_{z}$ and derive its fast GFT in Fig. 12 based on the centrosymmetry of the grid, characterized by the involution $\phi(i)=N+1-i$. Note that, if we flip the nodes $v_{9}$ to $v_{16}$ up to down, then ${\cal G}_{z}$ becomes a left-right symmetric grid, based on which we can derive ${\cal G}_{z}^{+}$ and ${\cal G}_{z}^{-}$ as in Fig. 12(a). The derived fast GFT diagram in Fig. 12(b) can thus provide a computational speedup for the coding framework in [47].

V-E Skeletal Graphs

In human action analysis, human body can be represented by a hierarchy of joints that are connected with bones. Human motion data can be obtained by two different types of techniques. First, motion capture is a process of directly extracting human movements from wearable devices, such as reflective markers attached near each joint. Second, using current methods such as OpenPose [48], skeletons can be obtained from images or videos in real time. With both techniques, the output of the system contains 3D coordinates of human joints. Then, we can consider the human skeleton as a graph, and motion data on the skeleton as graph signals. Further signal processing techniques can be applied to those signals to perform tasks such as classification and segmentation.

The work [49] has demonstrated that the GFT basis of the skeletal graph has localization properties useful in characterizing human motion. For example, the second GFT basis function has positive entries on joints in the upper body, and negative entries on those in the lower body. Thus, the resulting GFT coefficients can provide a discriminating power between different human actions.

Typical skeletal graphs are symmetric by construction, so a fast GFT can be obtained, as shown in Fig. 13, where a 15-node skeleton is considered. Such a fast GFT on skeletal graph can speed up the feature extraction procedure for further action classification tasks. Note that the butterfly stage is also available for skeletal graphs with non-uniform weights or different topologies, as long as the desired symmetry properties hold.

V-F Search of Symmetries in General Graphs

In the previous examples, symmetry properties of graphs can be easily identified by inspection. However, in general, and particularly for denser graphs, desired symmetry properties may not be straightforward to identify, or may not even exist. To design fast GFTs for graphs beyond the previous examples, an algorithm for searching a valid involution would thus be useful.

The number of involutions on $n$ elements is [36, Sec. 5.1.4]

[TABLE]

which asymptotically grows faster than a polynomial in $n$. This means that an exhaustive search of valid involutions among $T(n)$ possible candidates is a combinatorial problem. Here, we provide two methods to reduce the complexity of this search. More detailed illustrations and implementations of these methods can be found in [50].

V-F1 Pruning based on the degree list

Note that if ${\cal G}$ is $\phi$-symmetric, then the degrees of nodes $i$ and $\phi(i)$ must be equal for every $i$. This necessary condition for $\phi$-symmetry allows us to prune the involution search. In particular, we can compute the list of degrees first, then skip searching those involutions $\varphi$ with different degrees on nodes $i$ and $\varphi(i)$ for some $i$. For graphs with many distinct weight values, we tend to have many distinct node degrees, and thus the number of involutions that need to be searched can be significantly reduced.

V-F2 Searching of identical tree branches

Trees (i.e., graphs with no cycles) are connected graphs that have the smallest number of edges. This sparsity property implies that symmetry on trees can be characterized by pairs of identical subtrees (i.e., branches) whose roots are common or adjacent. For example, in Fig. 13, the two arms in the skeletal graph are identical branches that share a common root, and so are the two legs. Based on an algorithm proposed in [51], we provide in [50] an algorithm with ${\cal O}(n\log n)$ complexity that, for any given tree ${\cal G}$, finds all involutions $\phi$ such that ${\cal G}$ is $\phi$-symmetric.

V-G Complexity Analysis

In general, for a length-$n$ fast GFT with a layer of Haar units based on involution $\phi$, the number of multiplications is $p_{\phi}^{2}+(n-p_{\phi})^{2}$. This number is minimized to $n^{2}/2$ when $p_{\phi}=n/2$. This means that, in the best case scenario with one layer, the overall complexity is reduced by half, and the order of magnitude remains ${\cal O}(n^{2})$.

VI Experimental Results

In this section, we provide theoretical complexity analysis with the numbers of operations as well as experiments for empirical computation complexities of the fast GFTs. We have implemented several fast GFTs in C to simulate an environment closer to hardware. The source code for the experiments is available at [50].

VI-A Comparison with Matrix GFT

In the first experiment, we include the GFTs of several different graph topologies: the cycle graph with unit weights, the bi-diagonally symmetric 6-connected grid (as in Fig. 10, with $a=0.5$), the z-shaped grid with $w=2$, and the skeletal graph. For each graph we implement two fast GFTs with different sizes, and compare the runtime between the matrix GFT implementation and the fast GFT with butterfly stages. We include those GFTs in Figs. 9 to 13, together with larger graphs with the same topology types: the 80-node cycle graph, the 8$\times$8 bi-diagonally symmetric 4-connected grid, the 8$\times$8 z-shaped grid, and the 25-node skeletal graph used in [52]. Detailed design of their fast GFTs can be extended from the examples in Figs. 9 to 13. For each GFT, we generate 20000 graph signals with a proper length, whose entries are i.i.d. uniform random variables with range $[0,1]$. Then, we compute the percentage of runtime reduction for the symmetry-based fast GFT compared to the GFT realized by a single $n\times n$ matrix multiplication.

In Table III, we show, for each GFT, the numbers of additions (including subtractions), multiplications, and the empirical computation time reduction rate compared to matrix GFT in C implementation. We see that the fast GFT on skeletal graph in Fig. 13 with one butterfly stage leads to 45.5% speed improvement, and that on z-shaped grid in Fig. 12 gives around 41.5% runtime saving. Fast GFTs on cycle graphs and 6-connected grids that have multiple butterfly stages yield higher runtime reduction rates. From those results in Table III, we can see that the butterfly stages obtained from our proposed method lead to a significant speedup, and can be useful if the transform is required to be performed many times, and in a low-level or hardware implementation.

VI-B Comparison with Approximate Fast GFTs

In the second experiment, we compare our proposed method with an existing fast GFT approach [21] on graphs with symmetry properties. We consider two graphs for this experiment: the 8$\times$8 bi-diagonally symmetric grid with $a=0.5$, and the 8$\times$8 z-shaped grid with $w=2$. Note that, when the desired symmetry property is available, existing methods can be incorporated into the symmetry-based fast GFT scheme to speed up the computation of sub-GFTs such as ${\bf U}^{+}$ and ${\bf U}^{-}$. Thus, we can compare the following four GFT implementations:

Matrix GFT: an $n\times n$ matrix multiplication. 2. 2.

Haar + matrix GFT: symmetry-based fast GFT using Haar units, as shown in Figs. 10(b) and 12(b), where the sub-GFTs are implemented by full matrix multiplications. 3. 3.

Approximate GFT [21]: fast GFT using layers of Givens rotations found by the parallel truncated Jacobi algorithm–a greedy-based algorithm that progressively approximate $\hat{{\bf U}}^{\top}{\bf L}\hat{{\bf U}}$ to a diagonal matrix. The resulting GFT can be implemented using the schematic diagram as in Fig. 1. 4. 4.

Haar + approximate GFT: symmetry-based fast GFT with sub-GFTs implemented by approximate GFTs.

Let the GFT matrix associated to a GFT implementation be $\hat{{\bf U}}$, and the true GFT matrix be ${\bf U}$. We define two error metrics as follows.

Sign-normalized relative error (RE): we consider the relative error between two $n\times n$ orthogonal matrices,

[TABLE]

Note that if $\hat{{\bf U}}=-{\bf U}$, the RE will be large although they share a common eigen-structure. To avoid this sign ambiguity, we modify (23) by taking absolute values elementwise on $\hat{{\bf U}}^{\top}{\bf U}$:

[TABLE] 2. 2.

Empirical average error: let ${\cal X}=\{{\bf x}_{1},\dots,{\bf x}_{M}\}$ be the set of input signals, we define

[TABLE]

where ${\bf u}_{j}^{\top}{\bf x}_{i}$ and $\hat{{\bf u}}_{j}^{\top}{\bf x}_{i}$ are the $j$-th true and approximate GFT coefficients of ${\bf x}_{i}$. The absolute values are used to avoid the sign ambiguity.

For both graphs considered in this experiment, the eigenvalues of the Laplacians are all distinct, so the GFT bases have no rotation ambiguity.

We apply the method in [21] to obtain the parameters (angle and node pairings for Givens rotations) approximate GFTs, then implement the resulting fast algorithms in C, with different numbers of layers $J\in\{0,5,10,\dots\}$. We use $M=20000$ random samples as in Sec. VI-A, and obtain the error metrics, $\delta$ and $\epsilon$, for each GFT implementation.

The runtime versus sign-normalized RE, and versus empirical average error are shown in Figs. 14 and 15, respectively. We note that, first, the RE drops more steadily than the empirical error when the number of layers increases. This is related to the order of GFT basis functions. When more layers of Givens rotations are introduced, more GFT basis functions will be ordered correctly (i.e., smaller error in the final permutation operation \hbox{\boldmath\Pi}_{J+1} in Fig. 1). Indeed, we observe that when the number of correctly ordered GFT coefficients increases, the decrease of the empirical error is usually more significant than that of the relative error. The second observation is that for the two graphs with $n=64$ nodes, the approximate GFTs typically takes more than 20 layers to yield a sufficiently accurate GFT in terms of both error metrics. However, when more than 20 layers are used, the computation complexity becomes comparable or higher than Haar-matrix GFT, which provides exact GFT coefficients. Finally, we see that in both Figs. 14 and 15, the error of Haar + approximate GFT drops faster than that of approximate GFT. This means that by applying the symmetry property, we can obtain a significantly higher convergence rate for the approximation. This is a reasonable consequence, as our divide-and-conquer method reduces the dimension of the parameter estimation problem for [21].

VII Conclusion

In this paper, we have explored the relationship between the graph topology and properties in the corresponding GFT for fast GFT algorithm based on butterfly stages. We focus particularly on a component of the butterfly stage called Haar unit, and discuss the conditions for a stage of Haar units to be available in the GFT implementation. We have shown that a graph has a right butterfly stage with Haar units if it is k-regular bipartite. On the other hand, a left butterfly stage is available if the graph has certain symmetry properties. We have formally defined the relevant graph symmetry based on involution, i.e., pairing function of nodes. Then, we have proposed an approach, where once a graph symmetry is identified, we can decompose a graph ${\cal G}$ into two smaller graphs, ${\cal G}^{+}$ and ${\cal G}^{-}$, whose GFTs corresponds to the two parallel sub-transforms after a butterfly stage of Haar units. Again, from ${\cal G}^{+}$ and ${\cal G}^{-}$ we can explore subsequent butterfly stages if any desired symmetry property holds in them. Thus, this method enables us to explore butterflies stage by stage.

The desired symmetry properties typically arise in graphs that are nearly regular, symmetric by construction, or uniformly weighted. We have discussed several classes of those graphs: bipartite graphs, graphs with 2-sparse eigenvectors such as star and complete graphs, symmetric line and grid graphs, cycle graphs, and skeletal graphs. Relevant applications of those GFTs include video compression and human action analysis. Finally, we implement the fast GFT algorithms in C and compute the runtime saving for several graphs. The experiment results show that our method provides a significant computation time reduction compared to the GFT computed by matrix multiplication. It also outperforms existing fast approximate GFT approaches in terms of both complexity and accuracy for graphs with desired symmetry properties.

-A Proof of Theorem 4

In the following, we will repeatedly use (4) and (5). Also recall that in (17),

[TABLE]

From the block partition structure (12), we also have

[TABLE]

Such changes of indices will be used in the derivations below.

-A1 Edges of ${\cal G}^{+}$

With $i,j\in{\cal V}^{+}=\{1,\dots,n-p\}$ and $i\neq j$, we discuss three different cases separately, all based on (24). If $i,j\in{\cal V}_{X}=\{1,\dots,p\}$, then

[TABLE]

If $i\in{\cal V}_{X}=\{1,\dots,p\}$ and $j\in{\cal V}_{Z}=\{p+1,\dots,n-p\}$, then

[TABLE]

and the same holds for $i\in{\cal V}_{Z}$ and $j\in{\cal V}_{X}$. If $i,j\in{\cal V}_{Z}=\{p+1,\dots,n-p\}$, then

[TABLE]

-A2 Self-loops of ${\cal G}^{+}$

To express $s_{i}^{+}$, we discuss the cases with $i\in{\cal V}_{X}$ and $i\in{\cal V}_{Z}$ separately. If $i\in{\cal V}_{X}=\{1,\dots,p\}$, then from (24) we have

[TABLE]

If $i\in{\cal V}_{Z}=\{p+1,\dots,n-p\}$, then (24) gives

[TABLE]

-A3 Edges of ${\cal G}^{-}$

With $i,j\in{\cal V}_{Y}=\{n-p+1,\dots,n\}$ and $i\neq j$, from (25) we have

[TABLE]

-A4 Self-loops of ${\cal G}^{-}$

Here, we have $i\in{\cal V}_{Y}=\{n-p+1,\dots,n\}$, so, by (25),

[TABLE]

The results derived above apply to the case when $\phi=(n,n-1,\dots,1)$. For any arbitrary $\phi$, we can simply modify the sub-indices based on $\phi$. Thus, the results above can be written as in Lemma 4.