hide

# Natural Graph Wavelet Packet Dictionaries

Alexander Cloninger111Department of Mathematics and Halicioğlu Data Science
Institute, University of California, San Diego Haotian Li222Department of
Mathematics, University of California, Davis Naoki Saito††footnotemark:

###### Abstract

We introduce a set of novel multiscale basis transforms for signals on graphs
that utilize their “dual” domains by incorporating the “natural” distances
between graph Laplacian eigenvectors, rather than simply using the eigenvalue
ordering. These basis dictionaries can be seen as generalizations of the
classical Shannon wavelet packet dictionary and of classical time-frequency
adapted wavelet packets to arbitrary graphs, and does not rely on the
frequency interpretation of Laplacian eigenvalues. We describe the algorithms
(involving either vector rotations or orthogonalizations) to efficiently
approximate and compress signals through the best-basis algorithm, and
demonstrate the strengths of these basis dictionaries for graph signals on
sunflower graphs and road traffic networks.

##  1 Introduction and Motivation

There is an explosion of interest and demand to analyze data sampled on graphs
and networks. This has motivated development of more flexible yet
mathematically sound _dictionaries_ (i.e., an overcomplete collection of atoms
or basis vectors) for data analysis and signal processing on graphs. Our main
goal here is to build _smooth_ multiscale localized basis dictionaries on an
input graph, with beneficial reconstruction and sparsity properties, and to
fill the “gap” left from our previous graph basis dictionary constructions
[15, 14, 16, 17, 36] as we explain below. Our approach differs from the
standard literature as we fully utilize both the similarities between the
nodes (through the graph adjacency matrix) and the similarities between the
eigenvectors of the graph Laplacian matrix (through new non-trivial
eigenvector distances).

Previous approaches to construct such graph basis dictionaries break down into
two main categories. The first category partitions the nodes through recursive
graph cuts to generate multiscale basis dictionaries. This includes: the
Hierarchical Graph Laplacian Eigen Transform (HGLET) [15]; the Generalized
Haar-Walsh Transform (GHWT) [14]; its extension, the eGHWT [36]; and other
Haar-like graph wavelets (see, e.g., [29, 21, 10, 3, 38]). But their basis
vectors either are nonsmooth piecewise constants or have non-overlapping
supports. The second category uses spectral filters on the Laplacian (or
diffusion kernel) eigenvalues to generate multiscale smooth dictionaries. This
includes: the Spectral Graph Wavelet Transform [12]; Diffusion Wavelets [4];
extensions to spectral graph convolutional networks [22]. However, these
dictionaries ignore the fact that eigenvectors have a more complex
relationship than simply the difference in size of the corresponding
eigenvalues [34, 2, 24]. These relationships can result from eigenvector
localization in different clusters, differing scales in multi-dimensional
data, etc. These notions of similarity and difference between eigenvectors,
while studied in the eigenfunction literature [34, 2, 24], have yet to be
incorporated into building localized dictionaries on graphs.

We combine the benefits of both approaches to construct the graph equivalent
of spatial-spectral filtering. We have two approaches: one is to utilize the
_dual geometry_ of an input graph without partitioning the input graph, and
the other is to fully utilize clustering and partition information in both the
input graph and its dual domain.

Our first approach, detailed in Section 3, fills the “gap” in the cycle of our
development of the graph basis dictionaries, i.e., HGLET, GHWT, and eGHWT.
This approach is a direct generalization of the classical wavelet packet
dictionary [26, Chap. 8] to the graph setting: we hierarchically partition the
dual domain to generate a tree-structured “subbands” each of which is an
appropriate subset of the graph Laplacian eigenvectors. We also want to note
the following correspondence: The HGLET [15] is a graph version of the
Hierarchical Block DCT dictionary [26, Sec. 8.3] (i.e., the non-smooth non-
overlapping version of the local cosine dictionary [7], [26, Sec. 8.5]), and
the former exactly reduces to the latter if the input graph is $P_{N}$, a path
graph with $N$ nodes. The former hierarchically partitions the input graph
while the latter does the same (with a non-adaptive manner) on the unit
interval $[0,1]$ in the _time_ domain. On the other hand, the GHWT [14] is a
graph version of the Haar-Walsh wavelet packet dictionary [6], [26, Sec. 8.1],
and the former exactly reduces to the latter if the input graph is $P_{N}$.
The latter hierarchically partitions the interval $[0,N)$ in the _sequency_
domain while the former does the same by the graph domain partitioning plus
reordering; see [14, 16, 17] for the details. Our graph basis dictionary using
this first approach is a graph version of the _Shannon_ wavelet packet
dictionary [26, Sec. 8.1.2], which hierarchically partitions the interval
$[0,1/2)$ (or $[0,\pi]$ depending on how one defines the Fourier transform) in
the _frequency_ domain. Again, the former essentially reduces to the latter if
the input graph is $P_{N}$.

Our second approach, detailed in Section 4, partitions _both_ the input graph
_and_ its dual domain; more precisely, we first hierarchically partition the
dual domain, and then partition the input graph with constraints imposed by
the dual domain partition. This approach parallels and generalizes classical
time-frequency analysis, where the time domain is replaced by a general _node-
domain_ geometry and the frequency domain is replaced by a general
_eigenvector-domain_ organization that changes as a function of location in
the node-domain. A version of this approach of node-eigenvector organization
that maps the eigenvectors to a 1D projection of the graph has also been
considered as a visualization technique for low-frequency eigenfunctions on
clustered graphs [11].

We aim for the significance and impact of this research to be two fold. First,
these results will provide the first set of graph wavelet packet bases that
adaptively scale to the local structure of the graph. This is especially
important for graphs with complicated multiscale structure, whose graph
Laplacians have localized eigenfunctions, for example. This is an impactful
direction, as most of graph wavelet packet bases previously proposed only tile
the node-eigenfunction “plane” along the node “axis,” while Laplacian
eigenfunctions only tile that plane along the eigenfunction “axis”. Our
approach in Section 4 constructs filters in both the nodal domain and
eigenfunction domain, similar to the classical _time-frequency adapted_
wavelet packets that tile both the time and the frequency domains [39, 13].

Second, in the long term, this is a first method of systematically using the
novel concept of eigenvector dual geometry [34, 2, 24]. This direction can set
a path for future modification of spectral graph theory applications to
incorporate dual geometry.

Alex: I don’t think we need to add the following paragraph. Because we haven’t
talked about frames at all up to this point, and the PC-NGWP in Section 4 also
generates the basis dictionary, and needs to use the best-basis algorithm to
select an appropriate ONB. The graph wavelet basis is a particular choice from
this dictionary; and the same can be said for the VM-NGWP dictionary. _Earlier
Alex said: Check for rewording or moving to another place:_ The focus of this
paper is on construction of basis dictionaries, rather than overcomplete
frames. In Section 3 this is done through a best-basis algorithm for a given
signal, and in Section 4 this is done through construction of a signal
independent orthonormal wavelet basis. For this reason, in Section 5 we
compare to other basis dictionaries for a variety of signals.

The structure of this article is organized as follows. Section 2 reviews
fundamentals: the basics of graphs, i.e., graph Laplacians and graph Fourier
transform as well as graph wavelet transforms and frames that were proposed
previously. It also reviews non-trivial metrics of graph Laplacian
eigenvectors, which are used to construct the dual domain of an input graph.
Section 3 presents a natural graph wavelet basis constructed through
hierarchical partition of the dual graph of Laplacian eigenvectors. Section 4
presents a second version of a natural graph wavelet packet dictionary
constructed through a pair of hierarchical partitions, one on the input graph
and one on its dual domain. In Section 5, we demonstrate the usefulness of our
propose graph wavelet packet dictionaries in graph signal approximation using
numerical experiments. Code scripts can be found at [23]. Finally, we discuss
our findings gained through these numerical experiments and near-future
projects for further improvements of our dictionaries.

##  2 Background

###  2.1 Graph Laplacians and Graph Fourier Transform

Let $G=G(V,E)$ be an undirected connected graph. Let
${V=V(G)=\\{v_{1},v_{2},\ldots,v_{N}\\}}$ denote the set of nodes (or
vertices) of the graph, where $N\,:=\,|V(G)|$. For simplicity, we typically
associate each vertex with its index and write $i$ in place of $v_{i}$.
$E=E(G)=\\{e_{1},e_{2},\ldots,e_{M}\\}$ is the set of edges, where each
$e_{k}$ connects two vertices, say, $i$ and $j$, and $M\,:=\,|E(G)|$. In this
article we consider only finite graphs (i.e., $M,N<\infty$). Moreover, we
restrict to the case of simple graphs; that is, graphs without loops (an edge
connecting a vertex to itself) and multiple edges (more than one edge
connecting a pair of vertices). We use $\bm{f}\in\mathbb{R}^{N}$ to denote a
graph signal on $G$, and we define
$\bm{1}\,:=\,(1,\ldots,1)^{\scriptscriptstyle{\mathsf{T}}}\in\mathbb{R}^{N}$.

We now discuss several matrices associated with graphs. The information in
both $V$ and $E$ is captured by the _edge weight matrix_
$W(G)\in\mathbb{R}^{N\times N}$, where $W_{ij}\geq 0$ is the edge weight
between nodes $i$ and $j$. In an unweighted graph, this is restricted to be
either $0$ or $1$, depending on whether nodes $i$ and $j$ are adjacent, and we
may refer to $W(G)$ as an _adjacency matrix_. In a weighted graph, $W_{ij}$
indicates the _affinity_ between nodes $i$ and $j$. In either case, since $G$
is undirected, $W(G)$ is a symmetric matrix. We then define the _degree
matrix_ $D(G)$ as the diagonal matrix with entries $D_{ii}=\sum_{j}W_{ij}$.
With this in place, we are now able to define the _(unnormalized) Laplacian_
matrix, _random-walk normalized Laplacian_ matrix, and _symmetric normalized
Laplacian_ matrix, respectively, as

$\displaystyle L(G)$ $\displaystyle\,:=\,D(G)-W(G)$ $\displaystyle
L_{\mathrm{rw}}(G)$ $\displaystyle\,:=\,D(G)^{-1}L(G)$ (1) $\displaystyle
L_{\mathrm{sym}}(G)$ $\displaystyle\,:=\,D(G)^{-1/2}L(G)D(G)^{-1/2}.$

We use $0=\lambda_{0}\leq\lambda_{1}\leq\ldots\leq\lambda_{N-1}$ to denote the
sorted Laplacian eigenvalues and
$\bm{\phi}_{0},\bm{\phi}_{1},\ldots,\bm{\phi}_{N-1}$ to denote their
corresponding eigenvectors, where the specific Laplacian matrix to which they
refer will be clear from either context or subscripts. Denoting
$\Phi\,:=\,[\bm{\phi}_{0},\ldots,\bm{\phi}_{N-1}]$ and
$\Lambda\,:=\,\mathrm{diag}([\lambda_{0},\ldots,\lambda_{N-1}])$, the
eigendecomposition of $L(G)$ can be written as
$L(G)=\Phi\Lambda\Phi^{\scriptscriptstyle{\mathsf{T}}}$. Since we only
consider connected graphs here, we have $0=\lambda_{0}\lneqq\lambda_{1}$. This
second smallest eigenvalue $\lambda_{1}$ is called the _algebraic
connectivity_ of $G$ and the corresponding eigenvector $\bm{\phi}_{1}$ is
called the _Fiedler vector_ of $G$. The Fiedler vector plays an important role
in graph partitioning and spectral clustering; see, e.g., [41], which suggests
the use of the Fiedler vector of $L_{\mathrm{rw}}(G)$ for spectral clustering
over that of the other Laplacian matrices.

######  Remark 2.1.

In this article, we use the Fiedler vectors of $L_{\mathrm{rw}}$ of an input
graph and its subgraphs as a tool to hierarchically bipartition the graph
although any other graph partition methods can be used in our proposed
algorithms. However, note that we use the unnormalized graph Laplacian
eigenvectors of $L(G)$ for simplicity to construct the dual domain of $G$ and
consequently our graph wavelet packet dictionaries. In other words,
$L_{\mathrm{rw}}$ is only used to compute its Fiedler vector for our graph
partitioning purposes.

The graph Laplacian eigenvectors are often viewed as generalized Fourier modes
on graphs. Therefore, for any graph signal $\bm{f}\in\mathbb{R}^{N}$ and
coefficient vector $\bm{g}\in\mathbb{R}^{N}$, the _graph Fourier transform_
and _inverse graph Fourier transform_ [37] are defined by

$\displaystyle\mathcal{F}_{G}(\bm{f})\,:=\,\Phi^{\scriptscriptstyle{\mathsf{T}}}\cdot\bm{f}\in\mathbb{R}^{N}\quad\text{
and
}\quad\mathcal{F}^{\scriptscriptstyle{-1}}_{G}(\bm{g})\,:=\,\Phi\cdot\bm{g}\in\mathbb{R}^{N}.$
(2)

Since $\Phi$ is an orthogonal matrix, it is not hard to see that
$\mathcal{F}^{\scriptscriptstyle{-1}}_{G}\circ\mathcal{F}_{G}=I_{N}$. Thus, we
can use $\mathcal{F}_{G}$ as an analysis operator and
$\mathcal{F}^{\scriptscriptstyle{-1}}_{G}$ as a synthesis operator for graph
harmonic analysis.

###  2.2 Graph Wavelet Transforms and Frames

We now briefly review graph wavelet transforms and frames; see, e.g., [37, 31]
for more information. Translation and dilation are two important operators for
classical wavelet construction. However, unlike in $\mathbb{R}^{d}$
($d\in\mathbb{N}$) or its finite and discretized lattice graph
$P_{N_{1}}\times\cdots\times P_{N_{d}}$, we cannot assume the underlying graph
has self-symmetric structure in general (i.e., all the interior nodes do not
always have the same neighborhood structure, unless the graph is, for
instance, unweighted semiregular tree [1]). Therefore, it is difficult to
construct graph wavelet bases or frames by translating and dilating a single
mother wavelet function with a fixed shape, e.g., the Mexican hat mother
wavelet in $\mathbb{R}$, because the graph structure varies at different
locations. Instead, some researchers, e.g., Hammond et al. [12], constructed
wavelet frames by shifting smooth graph spectral filters to be centered at
different nodes. A general framework of building wavelet frames can then be
summarized as follows:

$\displaystyle\bm{\psi}_{j,n}\,:=\,\overbrace{\Phi
F_{j}\Phi^{\scriptscriptstyle{\mathsf{T}}}}^{\text{Filtering}}\bm{\delta}_{n}\quad\text{for
}j=0,1,\cdots,J\text{ and }n=1,2,\cdots,N,$ (3)

where the index $j$ stands for different scale of spectral filtering (the
greater $j$, the finer the scale), the index $n$ represents the center
location of the wavelet, the diagonal matrices $F_{j}\in\mathbb{R}^{N\times
N}$ satisfies $F_{0}(l,l)=h(\lambda_{l-1})$ and
$F_{j}(l,l)=g(s_{j}\lambda_{l-1})$ for $1\leq l\leq N$ and $1\leq j\leq J$.
Here, $h$ is a scaling function (which mainly deals with the small
eigenvalues), while $g$ is a graph wavelet generating kernel. For example, the
kernel proposed in [12] can be approximated by the Chebyshev polynomial and
lead to a fast algorithm. Note that $\\{s_{j}\\}_{1\leq j\leq J}$ are dilation
parameters and $\bm{\delta}_{n}$ is the standard basis vector centered at node
$n$.

Furthermore, one can show that as long as the generalized partition of unity

$A\cdot I_{N}\leq\sum_{j=0}^{N}F_{j}\leq B\cdot I_{N},\qquad 0<A\leq B$ (4)

holds, $\\{\bm{\psi}_{j,n}\\}_{0\leq j\leq J;1\leq n\leq N}$ forms a _graph
wavelet frame_ , which can be used to decompose and recover any given graph
signals [12].

###  2.3 Non-Trivial Eigenvector Distances

However, one important drawback of the above method is that the construction
of the spectral filters $F_{j}$ solely depends on the eigenvalue distribution
(except some flexibility in choosing the filter pair $(h,g)$, and the dilation
parameters $\\{s_{j}\\}$) and does _not_ reflect how the eigenvectors
_behave_. For simple graphs such as $P_{N}$ and $C_{N}$ (a cycle graph with
$N$ nodes), the graph Laplacian eigenvectors are global sinusoids whose
frequencies can be simply read off from the corresponding eigenvalues, as
discussed in [16, 34]. Hence, the usual _Littlewood-Paley wavelet theory_
(see, e.g., [9, Sec. 4.2], [18, Sec. 2.4]) applies for those simple graphs.
Unfortunately, the graph Laplacian eigenvectors of a general graph can behave
in a much more complicated manner than those of $P_{N}$ or $C_{N}$ (as
discussed in [35, 30, 34, 2, 24]). They cannot be described and ordered simply
by the corresponding eigenvalues. Therefore, it will be problematic to design
graph wavelet by using spectral filters built solely upon eigenvalues and we
need to find a way to distinguish eigenvector behaviors. To quantify the
difference of behaviors between the eigenvectors, however, we cannot use the
usual $\ell^{2}$-distances among them since they all have the same value
$\sqrt{2}$ due to their orthonormality.

As a remedy to these issues, we measure the “behavioral” difference between
the eigenvectors as discussed in [2, 34, 24], one of which is used in our
numerical experiments in Section 5. We briefly summarize this particular
behavioral distance below. See [2, 34, 24] for more information on the other
eigenvector metrics.

####  2.3.1 Hadamard (HAD) product affinity

The _Hadamard product affinity_ [2] defines a similarity between Laplacian
eigenvectors in terms of average local correlations. A benefit of this measure
is that it can be computed quickly. Let $\bm{\phi}_{i},\bm{\phi}_{j}$ be
eigenvectors of the graph Laplacian $L(G)$. It was shown in
[steinerberger2017spectral] that for all $x\in V$,

$\mathrm{e}^{-tL}(\bm{\phi}_{i}\circ\bm{\phi}_{j})(x)=\sum_{y\in
V}p_{t}(x,y)(\bm{\phi}_{i}(y)-\bm{\phi}_{i}(x))(\bm{\phi}_{j}(y)-\bm{\phi}_{j}(x)),$

where $\bm{\phi}_{i}\circ\bm{\phi}_{j}$ is the Hadamard (i.e., entrywise)
product of $\bm{\phi}_{i}$ and $\bm{\phi}_{j}$, $p_{t}(x,y)$ is the $(x,y)$-th
entry of $\mathrm{e}^{-tL}$, and $t$ satisfies

$\displaystyle\mathrm{e}^{-t\lambda_{i}}+\mathrm{e}^{-t\lambda_{j}}=1$ (5)

This implies that, rather than computing the local correlation between
eigenvectors for all $x\in V$, one can simply compute the Hadamard product
between the eigenvectors and apply the heat kernel matrix $\mathrm{e}^{-tL}$
at time $t$ satisfying (5). When
$\mathrm{e}^{-tL}(\bm{\phi}_{i}\circ\bm{\phi}_{j})$ is similar in magnitude to
$\bm{\phi}_{i}\circ\bm{\phi}_{j}$, this implies the local correlation between
$\bm{\phi}_{i}$ and $\bm{\phi}_{j}$ is large. This yields a natural similarity
measure between eigenvectors

$\alpha(\bm{\phi}_{i},\bm{\phi}_{j})\,:=\,\frac{\|\mathrm{e}^{-tL}(\bm{\phi}_{i}\circ\bm{\phi}_{j})\|_{2}}{\|\bm{\phi}_{i}\circ\bm{\phi}_{j}\|_{2}},$

which measures the fraction of $\bm{\phi}_{i}\circ\bm{\phi}_{j}$ that projects
into the graph Laplacian eigenvectors with smaller eigenvalues (i.e., more
global low-frequency eigenvectors). We define the _Hadamard product distance_
by (*** the logarithm the inverse of the Hadamard product affinity (5), i.e.,
We define the _Hadamard pseudometric_ by the logarithm of the inverse of

$d_{\mathrm{HAD}}(\bm{\phi}_{i},\bm{\phi}_{j})\,:=\,1/\alpha(\bm{\phi}_{i},\bm{\phi}_{j}).$

We note that this “distance” does not satisfy the metric axioms, yet it still
quantifies the behavioral difference between the eigenvectors as shown in [2].

Instead of the usual $\ell^{2}$-distance, we use the _absolute gradient_ of
each eigenvector as its feature vector describing its behavior. More
precisely, let $Q(G)\in\mathbb{R}^{N\times M}$ be the _incidence matrix_ of an
input graph $G(V,E,W)$ whose $k$th column indicates the head and tail of the
$k$th edge $e_{k}\in E$. However, we note that we need to orient each edge of
$G$ in an arbitrary manner to form a directed graph temporarily in order to
construct its incidence matrix. For example, suppose $e_{k}$ joins nodes $i$
and $j$, then we can set either
$(Q_{ik},Q_{jk})=(-\sqrt{W_{ij}},\sqrt{W_{ij}})$ or
$(\sqrt{W_{ij}},-\sqrt{W_{ij}})$. Of course, we set $Q_{lk}=0$ for $l\neq
i,j$. It is easy to see that
$Q(G)\,Q(G)^{\scriptscriptstyle{\mathsf{T}}}=L(G)$.

We now define the _DAG pseudometric_ between $\bm{\phi}_{i}$ and
$\bm{\phi}_{j}$ by

$d_{\mathrm{DAG}}(\bm{\phi}_{i},\bm{\phi}_{j})\,:=\,\||\nabla_{G}|\bm{\phi}_{i}-|\nabla_{G}|\bm{\phi}_{j}\|_{2}\quad\text{where
$|\nabla_{G}|\bm{\phi}\,:=\,\operatorname{abs}.(Q^{\scriptscriptstyle{\mathsf{T}}}\bm{\phi})\in\mathbb{R}^{M}_{\geq
0}$},$ (6)

where $\operatorname{abs}.(\cdot)$ applies the absolute value in the entrywise
manner to its argument. We note that this quantity is not a metric but a
_pseudometric_ because the identity of discernible of the axioms of metric is
not satisfied. In order to see the meaning of this quantity, let us analyze
its square as follows.

$\displaystyle d_{\mathrm{DAG}}(\bm{\phi}_{i},\bm{\phi}_{j})^{2}$
$\displaystyle=\left\langle{|\nabla_{G}|\bm{\phi}_{i}-|\nabla_{G}|\bm{\phi}_{j}},{|\nabla_{G}|\bm{\phi}_{i}-|\nabla_{G}|\bm{\phi}_{j}}\right\rangle_{E}$
$\displaystyle=\left\langle{|\nabla_{G}|\bm{\phi}_{i}},{|\nabla_{G}|\bm{\phi}_{i}}\right\rangle_{E}+\left\langle{|\nabla_{G}|\bm{\phi}_{j}},{|\nabla_{G}|\bm{\phi}_{j}}\right\rangle_{E}-2\left\langle{|\nabla_{G}|\bm{\phi}_{i}},{|\nabla_{G}|\bm{\phi}_{j}}\right\rangle_{E}$
$\displaystyle=\lambda_{i}+\lambda_{j}-\sum_{x\in V}\sum_{y\sim
x}|\bm{\phi}_{i}(x)-\bm{\phi}_{i}(y)|\cdot|\bm{\phi}_{j}(x)-\bm{\phi}_{j}(y)|\quad\text{thanks
to $QQ^{\scriptscriptstyle{\mathsf{T}}}=L$}$

where $\langle\cdot,\cdot\rangle_{E}$ is the inner product over edges. The
last term of the formula can be viewed as _a global average of absolute local
correlation_ between eigenvectors. In this sense, this quantity is related to
the _Hadamard-product affinity_ between eigenvectors proposed by Cloninger and
Steinerberger [2]. Note that the computational cost is $O(M)$ for each
$d_{\mathrm{DAG}}(\cdot,\cdot)$ evaluation provided that the eigenvectors have
already been computed.

###  2.4 Graph Wavelet Packets

Instead of building the graph wavelet packet dictionary by graph wavelet
frames as summarized in Section 2.2, one could also accomplish it by
generalizing the classical wavelet packets to graphs. The classical wavelet
packet decomposition (or dictionary construction) of a 1D discrete signal is
obtained by passing it through a full binary tree of filters (each node of the
tree represents either low-pass filtered or high-pass filtered versions of the
coefficients entering that node followed by the subsampling operation) to get
a set of binary-tree-structured coefficients [8], [26, Sec. 8.1]. This basis
dictionary for an input signal of length $N$ has up to $N(1+\log_{2}N)$ basis
vectors (hence clearly redundant), yet contains more than $1.5^{N}$ searchable
orthonormal bases (ONBs) [8, 39]. For the purpose of efficient signal
approximation, the _best-basis algorithm_ originally proposed by Coifman and
Wickerhauser [8] can find the most desirable ONB (and the expansion
coefficients of the input signal) for such a task among such an immense number
of ONBs. The best-basis algorithm requires a user-specified cost function,
e.g., the $\ell^{p}$-norm ($0<p\leq 1$) of the expansion coefficients for
sparse signal approximation, and the basis search starts at the bottom level
of the dictionary and proceeds upwards, comparing the cost of the coefficients
at the children nodes to the cost of the coefficients at their parents nodes.
This best-basis search procedure only costs $O(N)$ operations provided that
the expansion coefficients of the input signal have already been computed.

In order to generalize the classical wavelet packets to the graph setting,
however, there are two main difficulties: 1) the concept of the frequency
domain of a given graph is not well-defined; and 2) the relation between the
Laplacian eigenvectors and sample locations are much more subtle on general
graphs. For 1), we propose to construct a _dual graph_ $G^{\star}$ of the
input graph $G$ and view it as the natural spectral domain of $G$, and use any
graph partition method to hierarchically bipartition $G^{\star}$ instead of
building low and high pass filters like the classical case. This can be viewed
as the _generalized Littlewood-Paley theory_. For 2), we propose a node-
eigenvector organization algorithm called the _pair-clustering algorithm_ ,
which implicitly provides a downsampling process on graphs; see Section 4 for
the details.

##  3 Natural Graph Wavelet Packets using Varimax Rotations

Given a graph $G=G(V,E,W)$ with $|V|=N$ and the non-trivial distance $d$
between its eigenvectors (e.g., $d_{\mathrm{DAG}}$ of Eq. (6)), we build a
complete dual graph $G^{\star}=G^{\star}(V^{\star},E^{\star},W^{\star})$ by
viewing the eigenvectors as its nodes,
$V^{\star}=\\{\bm{\phi}_{0},\ldots,\bm{\phi}_{N-1}\\}$, and the non-trivial
affinity between eigenvector pairs as its edge weights,
$W^{\star}_{ij}=1/d(\bm{\phi}_{i-1},\bm{\phi}_{j-1})$, $i,j=1,2,\cdots,N$.
Using $G^{\star}$ for representing the graph spectral domain and studying
relations between the eigenvectors is clearly more _natural_ and effective
than simply using the eigenvalue magnitudes, as [2, 34, 24] hinted at. In this
section, we will propose one of our graph wavelet packet dictionary
constructions solely based on hierarchical bipartitioning of $G^{\star}$.
Basic Steps to generate such a graph wavelet packet dictionary for $G$ are
quite straightforward:

Step 1:

_Bipartition the dual graph $G^{\star}$ recursively_ via any method, e.g.,
spectral graph bipartition using the _Fiedler vectors_ ;

Step 2:

_Generate wavelet packet vectors_ using the eigenvectors belonging to each
subgraph of $G^{\star}$ that are _well localized on $G$_.

Note that Step 1 corresponds to bipartitioning the frequency band of an input
signal using the _characteristic functions_ in the classical setting. Hence,
our graph wavelet packet dictionary constructed as above can be viewed as a
graph version of the _Shannon_ wavelet packet dictionary [26, Sec. 8.1.2]. We
now describe the details of each step below.

###  3.1 Hierarchical Bipartitioning of $G^{\star}$

Let $V^{\star}=V^{\star(0)}_{0}\,:=\,[\bm{\phi}_{0},\ldots,\bm{\phi}_{n-1}]$
be the node set of the dual graph $G^{\star}$, which are simply the set of the
eigenvectors of the unnormalized graph Laplacian matrix $L(G)$. Suppose we get
the _hierarchical bipartition tree_ of $V^{\star(0)}_{0}$ as shown in Figure
1.

$V^{\star(0)}_{0}$$V^{\star(1)}_{0}$$V^{\star(2)}_{0}$ $V^{\star(2)}_{1}$
$V^{\star(1)}_{1}$$V^{\star(2)}_{2}$ $V^{\star(2)}_{3}$ Figure 1: The
hierarchical bipartition tree of the dual graph nodes $V^{\star}\equiv
V^{\star(0)}_{0}$, which corresponds to the frequency domain bipartitioning
used in the classical wavelet packet dictionary.

Hence, each $V^{\star(j)}_{k}$ contains an appropriate subset of the
eigenvectors of $L(G)$. Let $N^{j}_{k}$ be the number of such eigenvectors
belonging to $V^{\star(j)}_{k}$. With a slight abuse of notation, let
$V^{\star(j)}_{k}$ also denote a matrix of size $N\times N^{j}_{k}$ whose
columns are those eigenvectors. As we mentioned earlier, any graph
bipartitioning method can be used to generate this hierarchical bipartition
tree of $G^{\star}$. Typically, we use the Fiedler vector of the random-walk
normalized graph Laplacian matrix $L_{\mathrm{rw}}$ (see Eq. (2.1)) of each
subgraph of $G^{\star}$, whose use is preferred over that of $L$ or
$L_{\mathrm{sym}}$ as von Luxburg discussed [41].

###  3.2 Localization on $G$ via Varimax Rotation

For realizing Step 2 of the above basic algorithm, we propose to use the
_varimax rotation_ on $V^{\star(j)}_{k}\in\mathbb{R}^{N\times N^{j}_{k}}$ for
each $j$ and $k$. A varimax rotation is an orthogonal rotation, originally
proposed by Kaiser [20] and often used in _factor analysis_ (see, e.g., [28,
Chap. 11]), to maximize the variances of energy distribution (or a scaled
version of the _kurtosis_) of the input column vectors, which can also be
interpreted as the _approximate entropy minimization of the distribution of
the eigenvector components_ [33, Sec. 3.2]. For the implementation of the
varimax rotation algorithm, see Appendix A, which is based on the Basic
Singular Value (BSV) Varimax Algorithm of [19]. Thanks to the orthonormality
of columns of $V^{\star(j)}_{k}$, this is equivalent to finding an orthogonal
rotation that maximizes the overall _4th order moments_ , i.e.,

$\Psi^{(j)}_{k}\,:=\,V^{\star(j)}_{k}\cdot R^{(j)}_{k},\quad\text{where
$R^{(j)}_{k}=\arg\max_{R\in\mathrm{SO}(N^{j}_{k})}\sum_{p=1}^{N}\sum_{q=1}^{N^{j}_{k}}\left[\left(V^{\star(j)}_{k}\cdot
R\right)^{4}\right]_{p,q}$}.$ (7)

The column vectors of $\Psi^{(j)}_{k}$ are _more “localized” in the original
domain $G$_ than those of $V^{\star(j)}_{k}\equiv\Phi^{(j)}_{k}$. This type of
localization is important since the graph Laplacian eigenvectors in
$\Phi^{(j)}_{k}$ are of _global_ nature in general. We also note that the
column vectors of $\Psi^{(j)}_{k}$ are orthogonal to those of
$\Psi^{(j^{\prime})}_{k^{\prime}}$ as long as the latter is neither a direct
ancestor nor a direct descendant of the former. Hence, Steps 1 and 2 of the
above basic algorithm truly generate the graph wavelet packet dictionary for
an input graph signal, where one can run the best-basis algorithm of Coifman-
Wickerhauser [8] to extract the ONB most suitable for a task at hand (e.g., an
efficient graph signal approximation) once an appropriate cost function is
specified (e.g., the $\ell^{p}$-norm minimization, $0<p\leq 1$). Note also
that it is easy to extract a graph Shannon wavelet basis from this dictionary
by specifying the appropriate nodes in the binary tree structured dual graph
nodes explicitly, i.e.,
$\Psi^{(1)}_{1},\Psi^{(2)}_{1},\ldots,\Psi^{({j_{\mathrm{max}}})}_{1}$, and
the father wavelet vectors $\Psi^{({j_{\mathrm{max}}})}_{0}$. We refer to the
graph wavelet packet dictionary
$\left\\{\Psi^{(j)}_{k}\right\\}_{j=0:{j_{\mathrm{max}}};\,k=0:2^{j}-1}$
generated by this algorithm as the _Varimax Natural Graph Wavelet Packet_ (VM-
NGWP) dictionary.

Let us now demonstrate that our algorithm actually generates the classical
_Shannon_ wavelet packets dictionary [26, Sec. 8.1.2] when an input graph is
the simple path $P_{N}$.

(a) Father wavelet vectors $\Psi^{(4)}_{0}$

(b) Mother wavelet vectors $\Psi^{(4)}_{1}$

(c) Wavelet packet vectors $\Psi^{(4)}_{4}$

Figure 2: Some of the Shannon wavelet packet vectors on $P_{512}$

Note that the varimax rotation algorithm does not necessarily sort the vectors
as shown in Figure 2 because the minimization in Eq. (7) is the same modulo to
any permutation of the columns and any sign flip of each column. In other
words, to produce Figure 2, we carefully applied sign flip to some of the
columns, and sorted the whole columns so that each subfigure simply shows
translations of the corresponding wavelet packet vectors.

###  3.3 Computational Complexity

The varimax rotation algorithm of Appendix A is of iterative nature and is an
example of the BSV algorithms [19]: for each iteration at the dual node set
$V^{\star(j)}_{k}$, it requires computing the full Singular Value
Decomposition (SVD) of a matrix of size $N^{j}_{k}\times N^{j}_{k}$
representing a gradient of the objective function, which itself is computed by
multiplying matrices of sizes $N^{j}_{k}\times N$ and $N\times N^{j}_{k}$. The
convergence is checked with respect to the relative error between the current
and previous gradient estimates measured in the nuclear norm (i.e., the sum of
the singular values). For our numerical experiments in Section 5, we set the
maximum iteration as 1000 and the error tolerance as $10^{-12}$. Therefore, to
generate $\Psi^{(j)}_{k}$ for each $(j,k)$, the computational cost in the
worst case scenario is
$O\left(c\cdot(N^{j}_{k})^{3}+N\cdot(N^{j}_{k})^{2}\right)$ where $c=1000$ and
the first term accounts for the SVD computation and the second does for the
matrix multiplication. For a perfectly balanced bipartition tree with
$N=2^{j_{\mathrm{max}}}$, we have $N^{j}_{k}=2^{{j_{\mathrm{max}}}-j}$,
$j=0,\ldots,{j_{\mathrm{max}}}$, $k=0,\ldots,2^{j}-1$. Hence we have:

$\sum_{k=0}^{2^{j}-1}(N^{j}_{k})^{2}=\sum_{k=0}^{2^{j}-1}2^{2({j_{\mathrm{max}}}-j)}=2^{2{j_{\mathrm{max}}}-2j}\cdot
2^{j}=N^{2}\cdot 2^{-j},$ (8)

and

$\sum_{k=0}^{2^{j}-1}(N^{j}_{k})^{3}=\sum_{k=0}^{2^{j}-1}2^{3({j_{\mathrm{max}}}-j)}=2^{3{j_{\mathrm{max}}}-3j}\cdot
2^{j}=N^{3}\cdot 2^{-2j}.$

Note that at the bottom level $j={j_{\mathrm{max}}}$, each node is a leaf
containing only one eigenvector, and there is no need to do any rotation
estimation and computation. Hence, the total worst case computational cost is,
summing the cost $O\left(c\cdot(N^{j}_{k})^{3}+N\cdot(N^{j}_{k})^{2}\right)$
from $j=0$ to ${j_{\mathrm{max}}}-1$, i.e., $O((2+4c/3)N^{3}-2N^{2}-4c/3N)$,
so after all, it is an $O(N^{3})$ algorithm. In practice, the convergence is
often achieved with less than 1000 iterations at each node except possibly for
the nodes with small $j$ whose $N^{j}_{k}$ is large. For example, when
computing the VM-NGWP dictionary for the path graph $P_{512}$ shown in Figure
2, the average number of iterations over all the dual graph nodes
$\left\\{V^{(j)}_{k}\right\\}_{j=0:{j_{\mathrm{max}}};\,k=0:2^{j}-1}$ was
$70.43$ with the standard deviation $107.95$.

##  4 Natural Graph Wavelet Packets using _Pair-Clustering_

Another way to construct a natural graph wavelet packet dictionary is to mimic
the convolution and subsampling strategy of the classical wavelet packet
dictionary construction: form a full binary tree of spectral filters (i.e.,
the dual graph $G^{\star}$) and then perform the filtering/downsampling
process based on the relations between the nodes and the eigenvectors of $G$.
In order to fully utilize such relations, we look for a coordinated pair of
partitions on $G$ and $G^{\star}$, which is realized by our _pair-clustering_
algorithm described below. We will first describe the one-level pair-
clustering algorithm and then proceed to the hierarchical version.

###  4.1 One-Level Pair-Clustering

Suppose we partition the dual graph $G^{\star}$ into $K>2$ clusters using any
method including the spectral clustering [41] as we used in the previous
section. Let $V^{\star}_{1},\ldots,V^{\star}_{K}$ be those mutually disjoint
$K$ clusters of the nodes $V^{\star}$ representing the graph Laplacian
eigenvectors, i.e., $\displaystyle
V^{\star}=\bigsqcup_{k=1}^{K}V^{\star}_{k}$, which is also often written as
$\displaystyle\bigoplus_{k=1}^{K}V^{\star}_{k}$. Denote the cardinality of
each cluster as $N_{k}\,:=\,|V^{\star}_{k}|$, $k=1,\ldots,K$, and we clearly
have $\displaystyle\sum_{k=1}^{K}N_{k}=N$. Then, we also partition the
original nodes $V$ into mutually disjoint $K$ clusters, $V_{1},\ldots,V_{K}$
with the constraint that $|V_{k}|=|V^{\star}_{k}|=N_{k}$, $k=1,\ldots,K$, and
the members of $V_{k}$ and $V^{\star}_{k}$ are as “closely related” as
possible. The purpose of partitioning $V$ is to select appropriate graph nodes
as sampling points around which the graph wavelet packet vectors using the
information on $V^{\star}_{k}$ are localized. With a slight abuse of notation,
let $V$ also represent a collection of the standard basis vectors in
$\mathbb{R}^{N}$, i.e., $V\,:=\,\\{\bm{\delta}_{1},\ldots,\bm{\delta}_{N}\\}$,
where $\bm{\delta}_{k}(k)=1$ and $0$ otherwise. In order to formalize this
constrained clustering of $V$, we define the _affinity measure_ $\alpha$
between $V_{k}$ and $V^{\star}_{k}$ as follows:

$\displaystyle\alpha(V_{k},V^{\star}_{k})\,:=\,\sum_{\bm{\delta}\in
V_{k},\bm{\phi}\in
V^{\star}_{k}}|\left\langle{\bm{\delta}},{\bm{\phi}}\right\rangle|^{2},$ (9)

where $\left\langle{\cdot},{\cdot}\right\rangle$ is the standard inner product
in $\mathbb{R}^{N}$. Note that
$\displaystyle\alpha(V,V^{\star})=\sum_{\bm{\delta}\in V,\bm{\phi}\in
V^{\star}}|\left\langle{\bm{\delta}},{\bm{\phi}}\right\rangle|^{2}=\sum_{\bm{\phi}\in
V^{\star}}\|\bm{\phi}\|^{2}=N$. Denote the feasible partition space as

$U(V;N_{1},N_{2},\ldots,N_{K})\,:=\,\left\\{(V_{1},V_{2},\ldots,V_{K})\,\biggl{|}\,\bigsqcup_{k=1}^{K}V_{k}=V;\,|V_{k}|=N_{k},k=1,\ldots,K\right\\}.$

Now we need to solve the following optimization problem for a given partition
of $\displaystyle V^{\star}=\bigsqcup_{k=1}^{K}V^{\star}_{k}$:

$(V_{1},\ldots,V_{K})=\arg\max_{U(V;N_{1},N_{2},\ldots,N_{K})}\sum_{k=1}^{K}\alpha(V_{k},V^{\star}_{k})$
(10)

This is a discrete optimization problem. In general, it is not easy to find
the global optimal solution except for the case $K=2$. For $K=2$, we can find
the desired partition of $V$ by the following greedy algorithm: 1) compute
$\alpha(\\{\bm{\delta}\\},V^{\star}_{1})-\alpha(\\{\bm{\delta}\\},V^{\star}_{2})$
for each $\bm{\delta}\in V$; 2) select the largest $N_{1}$ such
$\bm{\delta}$’s in $V$, set them as $V_{1}$, and set $V_{2}=V\setminus V_{1}$.

When $K>2$, we can find a _local optimum_ by the similar strategy as above: 1)
compute the values $\alpha(\\{\bm{\delta}\\},V^{\star}_{1})$ for each
$\bm{\delta}\in V$; 2) select the largest $N_{1}$ such $\bm{\delta}$’s, and
set them as $V_{1}$; 3) compute the values
$\alpha(\\{\bm{\delta}\\},V^{\star}_{2})$ for each $\bm{\delta}\in V\setminus
V_{1}$, select the largest $N_{2}$ such $\bm{\delta}$’s, and set them as
$V_{2}$; 4) repeat the above process to produce $V_{3},\ldots,V_{K}$. While
this greedy strategy does not reach the global optimum of Eq (10), we find
that empirically the algorithm attains a reasonably large value of the
objective function.

###  4.2 Hierarchical Pair-Clustering

In order to build a multiscale graph wavelet packet dictionary, we develop a
hierarchical (i.e., multilevel) version of the pair-clustering algorithm.
First, let us assume that the hierarchical bipartition tree of $V^{\star}$ is
already computed using the same algorithm discussed in Section 3.1. We now
begin with level $j=0$ where $V^{(0)}_{0}$ is simply
$V=\\{\bm{\delta}_{1},\bm{\delta}_{2},\cdots,\bm{\delta}_{N}\\}$ and
$V^{\star(0)}_{0}$ is
$V^{\star}=\\{\bm{\phi}_{0},\bm{\phi}_{1},\cdots,\bm{\phi}_{N-1}\\}$. Then, we
perform one-level pair-clustering algorithm ($K=2$) to get
$\left(V^{\star(1)}_{0},V^{(1)}_{0}\right)$ and then
$\left(V^{\star(1)}_{1},V^{(1)}_{1}\right)$. We iterate the above process to
generate paired clusters $\left(V^{\star(j)}_{k},V^{(j)}_{k}\right)$,
$j=0,\ldots,{j_{\mathrm{max}}}$, $k=0,\ldots,2^{j}-1$. Note that we force the
partition so that each of the terminal nodes,
$V^{\star({j_{\mathrm{max}}})}_{k}$, $V^{({j_{\mathrm{max}}})}_{k}$, contains
only a single entry. In other words, the hierarchical pair-clustering
algorithm eventually will pair each standard basis $\bm{\delta}_{l}$ with
certain eigenvector $\bm{\phi}_{m}$.

###  4.3 Generating the NGWP Dictionary

Once we generate two hierarchical bipartition trees
$\left\\{V^{(j)}_{k}\right\\}$ and $\left\\{V^{\star(j)}_{k}\right\\}$, we can
proceed to generate the NGWP vectors $\left\\{\Psi^{(j)}_{k}\right\\}$ that
are necessary to form an NGWP dictionary. For each $\bm{\delta}_{l}\in
V^{(j)}_{k}$, we first compute the orthogonal projection of $\bm{\delta}_{l}$
onto the span of $V^{\star(j)}_{k}$, i.e.,
$\operatorname{span}\left(\Phi^{(j)}_{k}\right)$ where $\Phi^{(j)}_{k}$ are
those eigenvectors of $L(G)$ belonging to $V^{\star(j)}_{k}$. Unfortunately,
$\Phi^{(j)}_{k}\left(\Phi^{(j)}_{k}\right)^{\scriptscriptstyle{\mathsf{T}}}\bm{\delta}_{l}$
and
$\Phi^{(j)}_{k}\left(\Phi^{(j)}_{k}\right)^{\scriptscriptstyle{\mathsf{T}}}\bm{\delta}_{l^{\prime}}$
are not mutually orthogonal for $\bm{\delta}_{l},\bm{\delta}_{l^{\prime}}\in
V^{(j)}_{k}$ in general. Hence, we need to perform orthogonalization of the
vectors
$\left\\{\Phi^{(j)}_{k}\left(\Phi^{(j)}_{k}\right)^{\scriptscriptstyle{\mathsf{T}}}\bm{\delta}_{l}\right\\}_{l}$.
We use the _modified Gram-Schmidt with $\ell^{p}(0<p<2)$ pivoting
orthogonalization_ (MGSLp) [5] to generate the orthonormal graph wavelet
packet vectors associated with $V^{\star(j)}_{k}$ (and hence also
$V^{(j)}_{k}$). This MGSLp algorithm listed in Appendix B tends to generate
localized orthonormal vectors because the $\ell^{p}$-norm pivoting promotes
sparsity. We refer to the graph wavelet packet dictionary
$\left\\{\Psi^{(j)}_{k}\right\\}_{j=0:{j_{\mathrm{max}}};\,k=0:2^{j}-1}$
generated by this algorithm as the _Pair-Clustering Natural Graph Wavelet
Packet_ (PC-NGWP) dictionary.

###  4.4 Computational Complexity

At each node $V^{\star(j)}_{k}$ of the hierarchical bipartition tree of the
dual graph $G^{\star}$, the orthogonal projection of the standard basis
vectors in $V^{(j)}_{k}$ onto $\operatorname{span}\left(\Phi^{(j)}_{k}\right)$
and the MGSLp procedure are the two main computational burden for our PC-NGWP
dictionary construction. The orthogonal projection costs
$O\left(N\cdot(N^{j}_{k})^{2}+N\cdot N^{j}_{k}\right)$ while the MGSLp costs
$O\left(2N\cdot(N^{j}_{k})^{2}\right)$. Hence, the dominating cost for this
procedure is $O\left(3N\cdot(N^{j}_{k})^{2}\right)$ for each $(j,k)$. And we
need to sum up this cost on all the tree nodes. Let us analyze the special
case of the perfectly balanced bipartition tree with $N=2^{j_{\mathrm{max}}}$
as we did for the VM-NGWP in Section 3.3. In this case, the bipartition tree
has $1+{j_{\mathrm{max}}}$ levels, and $N^{j}_{k}=2^{{j_{\mathrm{max}}}-j}$,
$k=0,\ldots,2^{j}-1$. So, for the $j$th level, using Eq. (8), we have
$O(3N^{3}\cdot 2^{-j})$. Finally, by summing this from $j=0$ to
${j_{\mathrm{max}}}-1$ (again, no computation is needed at the bottom level
leaves), the total cost for PC-NGWP dictionary construction in this ideal case
is: $O(3N^{3}\cdot(2-2/N))\approx O(6N^{3})$. So, it still requires $O(N^{3})$
operations; the difference from that of the VM-NGWP is the constants, i.e.,
$6$ (PC-NGWP) vs $1000\cdot 4/3$ (the worst case VM-NGWP).

##  5 Applications in Graph Signal Approximation

In this section, we demonstrate the usefulness of our proposed NGWP
dictionaries in efficient approximation of graph signals on two graphs, and
compare the performance with the other previously-proposed methods, i.e., the
global graph Laplacian eigenbasis; graph Haar basis; graph Walsh basis; GHWT
coarse-to-fine (c2f) best basis [14]; GHWT fine-to-coarse (f2c) best basis
[14]; and eGHWT best basis [36]. Note that the graph Haar basis is a
particular basis choosable from the GHWT f2c dictionary and the eGHWT
dictionary while the graph Walsh basis is choosable from both versions of the
GHWT dictionaries as well as the eGHWT; see [17, 36] for the details. We use
the $\ell^{1}$-norm minimization as the best-basis selection criterion for all
the best bases in our experiments. The edge weights of the dual graph
$G^{\star}$ are the reciprocals of the DAG pseudometric between the
corresponding eigenvectors as defined in Eq. (6). For a given graph $G$ and a
graph signal $\bm{f}$ defined on it, we decompose $\bm{f}$ into those
dictionaries first. Then, to measure the approximation performance, we sort
the expansion coefficients in the non-increasing order of their magnitude, and
use the top $k$ most significant terms to approximate $\bm{f}$ where $k$
starts from $0$ up to about $0.5\times N$, i.e., 50% of the total number of
terms. All of the approximation performance is measured by the relative
$\ell^{2}$ approximation error with respect to the fraction of coefficients
retained, which ranges from $0$ to $0.5$.

###  5.1 Images Sampled on the Sunflower Graph

We consider the so-called “sunflower” graph shown in Fig. 3. This particular
graph has $400$ nodes and each edge weight is set as the reciprocal of the
Euclidean distance between the endpoints of that edge. Consistently counting
the number of spirals in such a sunflower graph gives rise to the _Fibonacci
numbers_ , i.e., $0,1,1,2,3,5,8,13,21,34,55,\dots$. One may observe that 8
clockwise and 13 counterclockwise spirals are emanating from the center in
Fig. 3. See, e.g., [27, 40] and our code SunFlowerGraph.jl in [23], for
algorithms to construct such sunflower grids and graphs. We can also view such
a distribution of nodes as a simple model of the distribution of
_photoreceptors in mammalian visual systems_ due to cell generation and
growth; see, e.g., [32, Chap. 9].

Figure 3: Sunflower graph ($N=400$); node radii vary for visualization purpose
only

Such a viewpoint motivates us the following sampling scheme: we overlay the
sunflower graph on several parts of the standard Barbara image and sample the
grayscale value at each node using the bilinear interpolation. In particular,
we sample two different regions: her left eye and pants, where quite different
image features are represented, i.e., a piecewise-smooth image containing
oriented edges and a textured image with directional oscillatory patterns,
respectively.

First, let us discuss our approximation experiments on Barbara’s eye graph
signal, which are shown in Fig. 4.

(a) The sunflower graph overlaid on Barbara’s left eye

(b) Barbara’s left eye as an input graph signal

(c) Approximation performance of various methods

Figure 4: Barbara’s left eye region sampled on the sunflower graph nodes (a)
as a graph signal (b); the relative $\ell^{2}$ approximation errors by various
methods (c)

From Fig. 4(c), we observe the following: 1) the VM-NGWP best basis performed
best closely followed by the PC-NGWP best basis; 2) the global graph Laplacian
eigenbasis worked relatively well; and 3) those bases chosen from the
Haar/Walsh wavelet packet dictionaries did not perform well. Note that the
GHWT c2f best basis turned out to be the graph Walsh basis in this graph
signal.

These observations can be attributed to the fact that this Barbara’s eye graph
signal is not of piecewise-_constant_ nature: it is a piecewise-_smooth_ graph
signal. Hence, those methods using smoother dictionary vectors worked better.
Yet, smooth localized basis vectors in the NGWP dictionaries made a difference
in performance compared to the global graph Laplacian eigenbasis. In other
words, the localized basis vectors were important for efficiently capturing
the features of this graph signal.

In order to examine what kind of basis vectors were chosen as the best basis
to approximate this Barbara’s eye signal, we display the most important nine
VM-NGWP best basis vectors in Fig. 5. The corresponding PC-NGWP best basis
vectors are relatively similar; hence they are not shown here. We note that
most of these top basis vectors essentially work as oriented edge detectors
for the eye including the iris part while the basis vectors #2 and #7 take
care of shading and peripheral features of this region.

Figure 5: Nine most important VM-NGWP vectors (the DC vector not shown) for
Barbara’s eye

Now, let us discuss our second approximation experiments: Barbara’s pants
region as an input graph signal as shown in Fig. 6. The nature of this graph
signal is completely different from the eye region: it is dominated by
directional oscillatory patterns of her pants.

(a) The sunflower graph overlaid on Barbara’s pants/knee region

(b) Barbara’s pants region as an input graph signal

(c) Approximation performance of various methods

Figure 6: Barbara’s pants region sampled on the sunflower graph nodes (a) as a
graph signal (b); the relative $\ell^{2}$ approximation errors by various
methods (c)

From Fig. 6(c), we observe the following: 1) NGWP best bases and the eGHWT
best basis performed very well and their performance is quite close until
about $30$% of the coefficients retained. Then eventually, the eGHWT best
basis outperformed all the others; 2) the GHWT f2c best basis performed
relatively well behind those three bases in 1); 3) there is a substantial gap
in performance between those four bases and the rest, i.e., the graph Haar
basis, the graph Walsh basis (= the GHWT c2f best basis in this case too), and
the global graph Laplacian eigenbasis.

We knew that the eGHWT is known to be quite efficient in capturing oscillating
patterns as shown by Shao and Saito for the graph setting [36] and by Lindberg
and Villemoes for the classical non-graph setting [25]. Hence, it is a good
thing to observe that our NGWPs are competitive with the eGHWT for this type
of textured signal. On the other hand, the graph Haar basis, the graph Walsh
basis, and the global graph Laplacian eigenbasis do not have oriented
anisotropic basis vectors, and consequently they performed poorly.

As in the case of Barbara’s eye signal, we display the most important nine
NGWP best basis vectors for approximating Barbara’s pants signal in Fig. 7,
but this time we display the PC-NGWP best basis vectors since that best basis
performed slightly better than the VM-NGWP best basis vectors. The
corresponding VM-NGWP best basis vectors, however, are relatively similar. We
note that the majority of these basis vectors are of high-frequency nature
than those for the eye signal shown in Fig. 5, which reflect the oscillating
anisotropic patterns of her pants. The basis vector #1 takes care of shading
in this region while the basis vector #5 is localized around the lower right
peripheral region. The latter takes care of the strong spotty pattern (the
dark point surrounded by light gray points) around that region in the input
graph signal shown in Figure 6(b).

Figure 7: Nine important PC-NGWP vectors (the DC vector not shown) for
Barbara’s pants

###  5.2 Toronto Traffic Network

We obtained the road network data of the City of Toronto from its open data
portal333URL: https://open.toronto.ca/dataset/traffic-signal-vehicle-and-
pedestrian-volumes. Using the street names and intersection coordinates
included in the dataset, we construct the graph representing the road network
there with $N=2275$ nodes and $M=3381$ edges. Fig. 8(a) displays this graph.
As before, each edge weight was set as the reciprocal of the Euclidean
distance between the endpoints of that edge.

We analyze two graph signals on this street network. 1) spatial distribution
of the street intersections and 2) pedestrian volume measured at such
intersections. The graph signal 1) was constructed by counting the number of
the nodes within the disk of radius of about $5$ km centered at each node. In
other words, this is just a smooth version of histogram of the distribution of
street intersections computed with the overlapping circular bins of equal
size. The longest edge length measured in the Euclidean distance among all
these $3381$ edges was chosen as this radius of this disk, which is located at
the northeast corner of this graph as one can easily see in Figure 8(a). The
graph signal 2) is the most recent 8 peak-hour volume counts collected at
intersections (i.e., nodes in this graph) where there are traffic signals. The
data was typically collected between the hours of 7:30 am and 6:00 pm, over
the period of 03/22/2004–02/28/2018.

(a) A smooth spatial distribution of the street intersections

(b) Approximation performance of various methods

Figure 8: (a) A graph signal representing the smooth spatial distribution of
the street intersections on the Toronto street network (a). The horizontal and
vertical axes of this plot represent the longitude and latitude geo-
coordinates of this area, respectively. (b) The results of our approximation
experiments.

From Fig. 8(b), we observe that qualitative behaviors of these error curves
are similar to those of Barbara’s eye signal shown in Fig. 4(c). That is, 1)
NGWP best bases outperformed all the others and the difference between the VM-
NGWP and the PC-NGWP is negligible; 2) the global graph Laplacian eigenbasis
worked quite well; 3) those bases based on Haar-Walsh wavelet packet
dictionaries did not perform well.

In order to examine what kind of basis vectors were chosen to approximate this
smooth histogram of street intersections, we display the most important nine
VM-NGWP best basis vectors in Fig. 9. We note that these top basis vectors
exhibit different spatial scales: the basis vectors #4, #6, and #9 are rather
localized while the basis vectors #2, #3, and #7 are of medium scale trying to
characterize the sharp gradient of the density function toward the most
crowded lower middle region. Then, there are coarse scale basis vectors such
as #1 and #5. The basis vector #8 is interesting: it is more oscillatory than
the others, and tries to characterize the difference within the dense lower
middle region.

Figure 9: Nine important VM-NGWP vectors (the DC vector not shown) for street
intersection density data on the Toronto street map

Now, let us analyze the pedestrian volume data measured at the street
intersections as shown in Fig. 10(a), which is quite localized around the most
dense region in the lower middle section of the street graph. Fig. 10(b) shows
the approximation errors of various methods.

(a) Pedestrian volume graph signal

(b) Approximation performance of various methods

Figure 10: The pedestrian volume graph signal on the Toronto street network
(a); the results of our approximation experiments (b)

From Fig. 10(b), we observe the following: 1) the eGHWT best basis clearly
outperformed all the other methods; 2) the graph Haar basis and both versions
of the GHWT best bases performed relatively well closely followed by the VM-
NGWP best basis; 3) the PC-NGWP best basis was outperformed by the VM-NGWP
best basis and the other three bases; 4) the global bases such as the graph
Laplacian eigenbasis and the graph Walsh based did not perform well.

In order to examine the performance difference between the VM-NGWP best basis
and the PC-NGWP best basis, we display the most important nine basis vectors
in Figures 11 and 12. We note that the top VM-NGWP best basis vectors exhibit
different spatial scales: the basis vectors #2 and #9 are rather localized
while the basis vectors #1, #3, #4, #7, #8 are of medium scale trying to
characterize the pedestrian data around the most crowded lower middle region.
Then, there are coarse scale basis vectors such as #1 and #5. The basis
vectors #5 and #6 try to delineate the boundary between the high and the low
pedestrian volume regions.

Figure 11: Nine important VM-NGWP vectors (the DC vector not shown) for
pedestrian volume data on the Toronto street map

On the other hand, the top PC-NGWP best basis vectors are much more localized
than those of the VM-NGWP best basis vectors. As one can see from Fig. 12, the
basis vectors #5, #8, #9 are extremely localized, and there are neither medium
nor coarse scale basis vectors in these top 9 basis vectors. The reason behind
these performance difference between the VM-NGWP and the PC-NGWP is the
following. The VM-NGWP best basis for this graph signal turned out to be a
true graph Shannon wavelet basis with scale $j=4$, i.e., the subspaces
$V^{\star(4)}_{0}$, $V^{\star(4)}_{1}$, $V^{\star(3)}_{1}$,
$V^{\star(2)}_{1}$, $V^{\star(1)}_{1}$ were selected as its best basis. On the
other hand, the PC-NGWP best basis turned out to be a graph wavelet basis with
$j=1$, i.e., the subspaces $V^{\star(1)}_{0}$ and $V^{\star(1)}_{1}$. Hence,
the latter has only fine scale basis vectors. This is because the PC-NGWP
procedure promotes localization of its basis vectors too strictly with scale
$j=1$ by using the orthogonal projections of the standard basis vectors
supported on $V^{(1)}_{k}$ of $G$ onto the subspace $V^{\star(1)}_{k}$ of
$G^{\star}$. Since the pedestrian volume data is quite non-smooth and
localized, $\bm{\delta}$-like basis vectors with scale $j=1$ in PC-NGWP tend
to generate sparser coefficients, i.e., having a small number of large
magnitude coefficients with many negligible ones. Therefore, the best basis
algorithm with $\ell^{1}$-norm ends up favoring those $\bm{\delta}$-like basis
vectors in PC-NGWP.

Figure 12: Nine important PC-NGWP vectors (the DC vector not shown) for
pedestrian volume data on the Toronto street map

##  6 Discussion

In this article, we proposed two ways to construct graph wavelet packet
dictionaries that fully utilize the natural dual domain of an input graph: the
VM-NGWP and the PC-NGWP dictionaries. Then, using two different graph signals
on each of the two different graphs, we compared their performance in
approximating a given graph signal against our previous multiscale graph basis
dictionaries, such as the GHWT and eGHWT dictionaries, which include the graph
Haar and the graph Walsh bases. Our proposed dictionaries outperformed the
others on piecewise smooth graph signals, and performed reasonably well for
graph signals sampled on an image containing oriented anisotropic texture
patterns. On the other hand, our new dictionaries were beaten by the eGHWT on
the non-smooth and localized graph signal. One of the potential reasons for
such a behavior is the fact that our dictionaries are a direct generalization
of the “Shannon” wavelet packet dictionaries, i.e., their “frequency” domain
support is localized and well controlled while the “time” domain support is
not compact. In order to improve the performance of our dictionaries for such
non-smooth localized graph signals, we need to bipartition $G^{\star}$
recursively but _smoothly with overlaps_ , which may lead to a graph version
of the _Meyer wavelet packet dictionary_ [26, Sec. 7.2.2, 8.4.2], whose basis
vectors are more localized in the “time” domain than those of the Shannon
wavelet packet dictionary. In fact, it is interesting to investigate such
smooth partitioning with overlaps not only on $G^{\star}$ but also on $G$
itself since it may lead to the graph version of the local cosine basis
dictionary [7], [26, Sec. 8.4.3].

We also note the differences between the VM-NGWP and the PC-NGWP dictionaries:
the PC-NGWP may allow one to “pinpoint” a particular node where the basis
vectors should concentrate in a more explicit manner than the VM-NGWP does.
The latter requires one to examine the computed basis vectors if one wants to
know where in the input graph nodes they concentrate. On the other hand, the
difference in their computational costs is just a constant in $O(N^{3})$
operations. Hence, it is important to investigate how to reduce the
computational complexity in both cases. One such possibility is to use only
the first $N_{0}$ graph Laplacian eigenvectors with $N_{0}\ll N$. Clearly, one
cannot represent a given graph signal precisely with $N_{0}$ eigenvectors, but
this scenario may be acceptable for certain applications including graph
signal clustering, classification, and regression. Of course, it is of our
interest to investigate whether we can come up with faster versions of varimax
rotation algorithm and MGSLp algorithm, which forms one of our future research
projects.

Do you think we should state the following explicitly?

* Building graph spectrogram for visualizing the eigenvectors

I now think that we do not need to state the following:

* Generalize this with LDB for graph signal classification problems

* What are critical applications that NGWP dictionaries outperform the others?

## Acknowledgments

This research was partially supported by the US National Science Foundation
grants DMS-1819222, DMS-1912747, CCF-1934568, DMS-2012266; the US Office of
Naval Research grant N00014-20-1-2381; and Russell Sage Foundation Grant 2196.

##  Appendix A Varimax Rotation

The algorithm of varimax rotation is a basic singular value (BSV) algorithm as
Jennrich demonstrated in [19]. The algorithm in Eq. (7) can be summarized as
follows:

  1. (0)

Initialize an orthogonal rotation matrix $R$.

  2. (1)

Compute $\mathrm{d}f/\mathrm{d}R$, where $f$ is the objective function defined
in Eq. (7). Here $\mathrm{d}f/\mathrm{d}R$ is the matrix of partial
derivatives of $f$ at $R$.

  3. (2)

Find the singular value decomposition $U\Sigma V^{*}$ of
$\mathrm{d}f/\mathrm{d}R$.

  4. (3)

Replace $R$ by $UV^{*}$ and go to (1) or stop.

Jennrich also showed that under certain general conditions, this algorithm
converges to a stationary point from any initial estimate. The above algorithm
seems to be the standard varimax rotation algorithm available in many
packages, e.g., MATLAB ®444MATLAB is a registered trademark of The MathWorks,
Inc., R, etc.

Input: Full column rank input matrix whose columns to be rotated
$A\in\mathbb{R}^{N\times m}$ ($m\leq N$); maximum number of iteration steps
maxit (default value: $1000$); relative tolerance tol (default value: 1e-12)

Output: Rotated matrix $B\in\mathbb{R}^{N\times m}$

$B=A$

// initialize the output matrix

$S=0$

// initialize the nuclear norm

for  _$i$ from $1$ to maxit_ do

$S_{0}=S$

$[U,\Sigma,V]=\operatorname{svd}(A^{\scriptscriptstyle{\mathsf{T}}}\cdot(N\cdot
B^{\circ 3}-B\cdot\operatorname{diag}(B^{\scriptscriptstyle{\mathsf{T}}}\cdot
B)))$

// $B^{\circ 3}\,:=\,B\circ B\circ B$

// where $\circ$ is the Hadamard (entrywise) product

$T=U\cdot V^{*}$

// update the orthogonal rotation matrix

$S=\operatorname{trace}(\Sigma)$

$B=A\cdot T$

// update the rotated matrix

if  _$|S-S_{0}|/S <$ tol_ then

break

// stop when $S$ does not change much

end if

end for

return _$B$_

Algorithm 1 The Varimax Rotation Algorithm

##  Appendix B Modified Gram-Schmidt with $\ell^{p}$ Pivoting
Orthogonalization

We implemented a simplified version (i.e., Algorithm 2) of the modified Gram-
Schmidt with mixed $\ell^{2}$-$\ell^{p}$ ($0<p<2$) pivoting algorithm in [5].
Our version skips the the step of computing the largest $\ell^{2}$ norm and
picking the parameter $\lambda$ (a notation used in [5]) to increase the
numerical stability. Instead, we directly set up a tolerance parameter, i.e.,
$\mathrm{tol}$ in Algorithm 2, for the robustness. On the other hand, we keep
the $\ell^{p}$ ($0<p<2$) pivoting portion in MGS (i.e., always perform the
orthogonalization process of the vector with minimum $\ell^{p}$-norm in the
candidate pool), which nicely preserves the sparsity of the obtained wavelet-
like vectors after the orthogonalization process. The MGSLp algorithm is
summarized as follows.

Input: List of unit vectors $v=[v_{1},\ldots,v_{m}]\in\mathbb{R}^{N\times m}$;
norm parameter $0<p<2$ (default value: $1$); error tolerance tol (default
value: 1e-12)

Output: List of orthonormal vectors
$q=[q_{1},\ldots,q_{r}]\in\mathbb{R}^{N\times r}$ where
$r=\operatorname{rank}(v)$

$q=\emptyset$

// initialize the output list

$w=[\|v_{1}\|_{p},\ldots,\|v_{m}\|_{p}]$

for  _$i$ from $1$ to $m$_ do

$k=i-1+\operatorname{findmin}(w)$

// find the minimum $\ell^{p}$-norm index

swap($v_{i}$, $v_{k}$)

// pivoting

if  _$\|v_{i}\|_{2} <\textnormal{{tol}}$_ then

break

// check linear dependency

end if

$\tilde{v}=v_{i}/\|v_{i}\|_{2}$

$w=\emptyset$

// re-initialize the $\ell^{p}$-norm vector

for  _$j$ from $i+1$ to $m$_ do

$v_{j}=v_{j}-(\tilde{v}^{\scriptscriptstyle{\mathsf{T}}}\cdot v_{j})\tilde{v}$

$w$ = append($w$, $\|v_{j}\|_{p}$)

end for

$q$ = append($q$, $\tilde{v}$)

end for

return _$q$_

Algorithm 2 Modified Gram-Schmidt Orthogonalization with $\ell^{p}$ pivoting

## References

  * [1] T. Bıyıkoğlu, J. Leydold, and P. F. Stadler, Laplacian Eigenvectors of Graphs, vol. 1915, Springer, 2007.  Lecture Notes in Mathematics.
  * [2] A. Cloninger and S. Steinerberger, On the dual geometry of Laplacian eigenfunctions, Experimental Mathematics, (2018).  DOI:10.1080/10586458.2018.1538911.
  * [3] R. R. Coifman and M. Gavish, Harmonic analysis of digital data bases, in Wavelets and Multiscale Analysis: Theory and Applications, J. Cohen and A. I. Zayed, eds., Applied and Numerical Harmonic Analysis, Boston, MA, 2011, Birkhäuser, pp. 161–197.
  * [4] R. R. Coifman and M. Maggioni, Diffusion wavelets, Appl. Comput. Harmon. Anal., 21 (2006), pp. 53–94.
  * [5] , Diffusion wavelets, Appl. Comput. Harmon. Anal., 21 (2006), pp. 53–94.
  * [6] R. R. Coifman and Y. Meyer, Nouvelles bases orthonormées de $L^{2}(\mathbb{R})$ ayant la structure du systèm de Walsh, preprint, Dept. of Mathematics, Yale University, New Haven, CT, Aug. 1989.
  * [7] , Remarques sur l’analyse de Fourier à fenêtre, Comptes Rendus Acad. Sci. Paris, Série I, 312 (1991), pp. 259–261.
  * [8] R. R. Coifman and M. V. Wickerhauser, Entropy-based algorithms for best basis selection, IEEE Trans. Inform. Theory, 38 (1992), pp. 713–718.
  * [9] I. Daubechies, Ten Lectures on Wavelets, vol. 61 of CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM, Philadelphia, PA, USA, 1992.
  * [10] M. Gavish, B. Nadler, and R. R. Coifman, Multiscale wavelets on trees, graphs and high dimensional data: theory and applications to semi supervised learning, in Proc. 27th Intern. Conf. Machine Learning, J. Fürnkranz and T. Joachims, eds., Haifa, Israel, 2010, Omnipress, pp. 367–374.
  * [11] B. Girault and A. Ortega, What’s in a frequency: new tools for graph Fourier transform visualization, arXiv preprint arXiv:1903.08827, (2019).
  * [12] D. K. Hammond, P. Vandergheynst, and R. Gribonval, Wavelets on graphs via spectral graph theory, Appl. Comput. Harmon. Anal., 30 (2011), pp. 129–150.
  * [13] C. Herley, Z. Xiong, K. Ramchandran, and M. T. Orchard, Joint space-frequency segmentation using balanced wavelet packet trees for least-cost image representation, IEEE Trans. Image Process., 6 (1997), pp. 1213–1230.
  * [14] J. Irion and N. Saito, The generalized Haar-Walsh transform, in Proc. 2014 IEEE Workshop on Statistical Signal Processing, 2014, pp. 472–475.
  * [15] , Hierarchical graph Laplacian eigen transforms, JSIAM Letters, 6 (2014), pp. 21–24.
  * [16] , Applied and computational harmonic analysis on graphs and networks, in Wavelets and Sparsity XVI, Proc. SPIE 9597, M. Papadakis, V. K. Goyal, and D. Van De Ville, eds., 2015.  Paper # 95971F.
  * [17] , Efficient approximation and denoising of graph signals using the multiscale basis dictionaries, IEEE Trans. Signal and Inform. Process. Netw., 3 (2017), pp. 607–616.
  * [18] S. Jaffard, Y. Meyer, and R. D. Ryan, Wavelets: Tools for Science & Technology, SIAM, Philadelphia, PA, 2001.
  * [19] R. I. Jennrich, A simple general procedure for orthogonal rotation, Psychometrika, 66 (2001), pp. 289–306.
  * [20] H. F. Kaiser, The varimax criterion for analytic rotation in factor analysis, Psychometrika, 23 (1958), pp. 187–200.
  * [21] A. Lee, B. Nadler, and L. Wasserman, Treelets—an adaptive multi-scale basis for sparse unordered data, Ann. Appl. Stat., 2 (2008), pp. 435–471.
  * [22] R. Levie, F. Monti, X. Bresson, and M. M. Bronstein, CayleyNets: Graph convolutional neural networks with complex rational spectral filters, IEEE Trans. Signal Process., 67 (2018), pp. 97–109.
  * [23] H. Li, Natural graph wavelet dictionaries.  https://github.com/haotian127/Natural-Graph-Wavelet-Dictionaries, 2020.
  * [24] H. Li and N. Saito, Metrics of graph Laplacian eigenvectors, in Wavelets and Sparsity XVIII, Proc. SPIE 11138, D. Van De Ville, M. Papadakis, and Y. M. Lu, eds., 2019.  Paper #111381K.
  * [25] M. Lindberg and L. F. Villemoes, Image compression with adaptive Haar-Walsh tilings, in Wavelet Applications in Signal and Image Processing VIII, Proc. SPIE 4119, A. Aldroubi, A. F. Laine, and M. A. Unser, eds., 2000, pp. 911–921.
  * [26] S. Mallat, A Wavelet Tour of Signal Processing, Academic Press, Burlington, MA, third ed., 2009.
  * [27] A. M. Mathai and T. A. Davis, Constructing the sunflower head, Math. Biosci., 20 (1974), pp. 117–133.
  * [28] S. A. Mulaik, Foundations of Factor Analysis, CRC Press, 2nd ed., 2010\.
  * [29] F. Murtagh, The Haar wavelet transform of a dendrogram, J. Classification, 24 (2007), pp. 3–32.
  * [30] Y. Nakatsukasa, N. Saito, and E. Woei, Mysteries around the graph Laplacian eigenvalue $4$, Linear Algebra Appl., 438 (2013), pp. 3231–3246.
  * [31] A. Ortega, P. Frossard, J. Kovačević, J. M. F. Moura, and P. Vandergheynst, Graph signal processing: Overview, challenges, and applications, Proc. IEEE, 106 (2018), pp. 808–828.
  * [32] R. W. Rodieck, The First Steps in Seeing, Sinauer Associates, Inc., Sunderland, MA, 1998.
  * [33] N. Saito, The generalized spike process, sparsity, and statistical independence, in Modern Signal Processing, D. Rockmore and J. D. Healy, eds., vol. 46 of MSRI Publications, Cambridge University Press, 2004, pp. 317–340.
  * [34] , How can we naturally order and organize graph Laplacian eigenvectors?, in Proc. 2018 IEEE Workshop on Statistical Signal Processing, 2018, pp. 483–487.
  * [35] N. Saito and E. Woei, On the phase transition phenomenon of graph Laplacian eigenfunctions on trees, RIMS Kôkyûroku, 1743 (2011), pp. 77–90.
  * [36] Y. Shao and N. Saito, The extended generalized Haar-Walsh transform and applications, in Wavelets and Sparsity XVIII, Proc. SPIE 11138, D. Van De Ville, M. Papadakis, and Y. M. Lu, eds., 2019.  Paper #111380C.
  * [37] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, The emerging field of signal processing on graphs, IEEE Signal Processing Magazine, 30 (2013), pp. 83–98.
  * [38] A. D. Szlam, M. Maggioni, R. R. Coifman, and J. C. Bremer, Jr., Diffusion-driven multiscale analysis on manifolds and graphs: top-down and bottom-up constructions, in Wavelets XI, Proc. SPIE 5914, M. Papadakis, A. F. Laine, and M. A. Unser, eds., 2005.  Paper # 59141D.
  * [39] C. M. Thiele and L. F. Villemoes, A fast algorithm for adapted time-frequency tilings, Appl. Comput. Harmon. Anal., 3 (1996), pp. 91–99.
  * [40] H. Vogel, A better way to construct the sunflower head, Math. Biosci., 44 (1979), pp. 179–189.
  * [41] U. von Luxburg, A tutorial on spectral clustering, Stat. Comput., 17 (2007), pp. 395–416.


# VU-BERT: A Unified framework for Visual Dialog

###### Abstract

The visual dialog task attempts to train an agent to answer multi-turn
questions given an image, which requires the deep understanding of
interactions between the image and dialog history. Existing researches tend to
employ the modality-specific modules to model the interactions, which might be
troublesome to use. To fill in this gap, we propose a unified framework for
image-text joint embedding, named VU-BERT, and apply patch projection to
obtain vision embedding firstly in visual dialog tasks to simplify the model.
The model is trained over two tasks: masked language modeling and next
utterance retrieval. These tasks help in learning visual concepts, utterances
dependence, and the relationships between these two modalities. Finally, our
VU-BERT achieves competitive performance (0.7287 NDCG scores) on VisDial v1.0
Datasets.

Index Terms—  Multi-Modal, Visual Dialog, Patch Embedding, Transformer

##  1 Introduction

Interactions between humans and the world generally involve joint perception
of information in miscellaneous modalities, _e.g._ , image, text and audio. On
the contrary, most of the well-known machine learning models are modality-
specific like BERT [1, 2] or ResNet [3], agnostic to the information from
interweaving of different modalities. So motivated, multi-modal tasks have
recently gained increasing popularity, especially in the fields of vision and
language. At present, popular visual and language tasks include Visual Caption
(VC) [4, 5], Visual Grounding [6, 7], Visual Question Answering (VQA) [4, 7,
8] and Visual Dialog (VD) [9, 10, 11]. VQA attempts to predict a correct
answer to questions given some background texts and images. The task of VD is
firstly proposed by [12], where a dialog agent needs to answer a series of
questions grounded by an image. Compared to VQA, VD is more challenging as the
model needs not only to present high-level semantic understanding according to
image and dialog history, but also to process and reason through the longer
contexts than VQA.

To address VD tasks, existing works generally encode the input contexts and
images separately with modality-specific modules, then combine the information
from the two sources with a variety of cross-modal attention mechanisms, as in
DVAN [13], DualVD [14], VisDial-BERT [15], _etc_. Different from them, our VU-
BERT leverages a single-stream transformer-based model to learn interactions,
which jointly encodes information from the images and text without modality-
specific modules. Most closely to our VU-BERT is VD-BERT [11], which also
models interactions by a unified transformer. However, VD-BERT must
concatenate every answer candidate with the input and go through a forward
pass of the entire model, which requires high computing resources and cannot
scale to a large number of candidates like ours. Besides, our VU-BERT adapts
patch embedding firstly in VD tasks while all of prior models process images
complicatedly and inefficiently by ResNet [3] or Faster R-CNN [16].

To explore the interactions between images and dialog history, our framework
utilizes BERT [1] architecture, and trains the model over two objectives:
masked language modeling and next utterance retrieval. Different from previous
work, we adopt a tractable visual embedding introduced in ViT [17], which
yields significant reduction on latency time and parameter scale [4]. The key
contributions of our study can be summarized as follows:

  * •

We propose a unified framework VU-BERT for effective vision and dialog
representations to predict the next utterance by selecting or generating
responses. And we achieve competitive performance on visual dialog tasks.

  * •

Different from existing works [9, 10, 11], we attempt to obtain visual
embedding by patch projection to the VD tasks, which leads to a simple
architecture without any region features or deep convolution operations.

  * •

We conduct extensive experiments to verify the effectiveness of our framework.
According to the ablation studies, we found that different dialog turns affect
the performance, which indicates that it is not suitable to apply models
specific to VQA or VC tasks directly to VD tasks.

Fig. 1: Overview of our VU-BERT. The image $I$ is split to patches and
linearly projected to obtain patch embedding. The input consists of image and
text embedding via summing up its position embedding and segment embedding. We
keep the input as $(I,C,H,U_{m})$ and the response candidates $A$ respectively
passing through a transformer encoder.

##  2 methodology

In this section, we discuss our visual dialog model VU-BERT, including two
novel approaches aimed at capturing better representations of vision and
utterances as shown in Fig 1. We first employ a unified vision-dialog
Transformer to encode both the image and dialog history. Then, we adopt
visually grounded masked language modeling for effective interactions between
visual and textual features and next utterance retrieval for enhancing
utterances dependence to train the model. Next, we describe the components of
this model in detail.

###  2.1 Embeddding

Dialog Embedding. Given a $m$-th turn dialog represented by a list of
utterances $U=[u_{1},\ldots,u_{m}]$ grounded in the image $I$, an utterance is
first split into words $u_{i}=[w_{1},\ldots,w_{n}]$ with length of $n$ by the
same WordPiece tokenizer [18]. We follow BERT and each sub-word is obtained
via summing up its word embedding, position embedding and segment embedding.

Vision Embedding. Different from extracting object-level features as vision
features [9, 10, 11], we use ViT[17] as a backbone network to process images,
which is faster than object detectors. According to ViT, Patch Embedding will
split the height and width of the input image $I\in R^{H\times W\times C}$ to
$N=\frac{HW}{P^{2}}$ patches according P, and then flatten and reshape the
patches to $v\in R^{N\times\left(P^{2}\times C\right)}$ through linear
transformation. Similar to word embedding in BERT, we add segment embedding
and position embedding to patch embedding as visual embedding where all the
patches share the same segment embedding.

###  2.2 Network Architecture

VU-BERT, which unifies the multi-model architecture into a single model, is
described in Fig.1. We prepare the input sequence as $x=\\{[CLS]$, $v_{1}$,
$\ldots$, $v_{t}$, $[EOI],C$, $[SEP]$, $u_{1}$, $[SEP],\ldots,u_{m},[END]\\}$.
The special token $[CLS]$ and $[END]$ denote the beginning and the end of the
sequence respectively and $[SEP]$ separates the different utterances. In
addition, we add the token $[EOI]$ to indicate the end of image input.
Moreover, we combine the input sequence with its position embedding and
segment embedding. Then we encode the sequences into cross-modal
representations $H=[h_{cls},h_{1},\ldots,h_{end}]$ using the framework of
Transformer. According to the seq2seq self-attention masks [19, 20], we mask
out all next-turn utterance tokens by setting $-\infty$ in $M$ to get auto-
regressive attentions while keeping tokens in the image $I$, caption $C$,
dialog history $H$ and current utterance $U_{m}$ visible, which can
recursively generate an answer by predicting the masked tokens. Similarly, we
train VU-BERT with the NUR loss using bidirectional self-attention masks to
select a correct response from response candidates.

###  2.3 Training Objectives

To obtain a strong and general representation for the model, we train VU-BERT
with the following two losses:

Table 1: Performance of discriminative/generation models on VisDial v0.9 dataset. Higher is better for MRR and recall@k, while lower is better for Mean. Model | MRR $\uparrow$ | R@1 $\uparrow$ | R@5 $\uparrow$ | R@10 $\uparrow$ | Mean $\downarrow$
---|---|---|---|---|---
LF | 0.5807/0.5199 | 43.82/41.83 | 74.68/61.78 | 84.07/67.59 | 5.78/17.07
HRE | 0.5846/0.5237 | 44.67/42.29 | 74.50/62.18 | 84.22/67.92 | 5.72/17.07
DVAN | 0.6675/0.5594 | 53.62/46.58 | 82.85/65.50 | 90.72/71.25 | 3.93/14.79
FGA | 0.6525/0.4138 | 51.43/27.42 | 82.08/56.33 | 89.56/71.32 | 4.35/9.10
VD-BERT | 0.7004/0.5595 | 55.79/46.83 | 85.34/65.43 | 92.68/72.05 |  4.04/13.18
VU-BERT | 0.6333/0.5403 | 48.71/44.50 | 81.03/62.60 | 89.10/71.70 | 4.19/12.49

####  2.3.1 Masked Language Modeling (MLM)

Similar to MLM in BERT, each token in utterance is randomly masked out with
15% probability and replaced with a special token [MASK]. The model is trained
to predict the masked words, based on the unmasked words and the visual
features. The task drives the network to not only model the dependencies in
sentence words, but also to align the visual and utterance contents. We denote
the visual embedding as $v={v_{1},\ldots,v_{t}}$ and the mask indices as $k$,
where $w_{\backslash k}$ represents as $w_{\backslash
k}=\\{w_{1},\ldots,w_{k-1},[MASK],w_{k+1},\ldots,w_{n}\\}$. The parameters are
learned to minimize the cross-entroy loss computed using the predicted tokens
and the original tokens. Then the MLM loss function is shown in Eq (1):

$\displaystyle L_{mlm}=-E_{(w,v)\sim D}\log P(w_{k}|w_{\backslash k},v)$ (1)

####  2.3.2 Next Utterance Retrieval (NUR)

Given the utterances $[u_{1},\ldots,u_{m-1}]$, the aim of NUR is to select the
correct next utterance $u_{m}$ from a set of $N$ candidate responses. For this
task, we need to respectively encode entire dialog context and each candidate
response by using BERT, which allows for the precomputation of all possible
candidates and can scale to millions of candidates. A list of candidate
responses can be created automatically, with the positive sample being the
next utterance in the dialog history, and making the assumption that randomly
sampled utterances from different dialogues are likely to be negatives. We
calculate the probability of $u_{i}$ being the true next utterance and
minimize the cross-entropy loss of the next utterance, which is analogous to
language modeling. The NUR loss function is:

$\displaystyle L_{nur}$ $\displaystyle=-\log p(u_{i}|h_{cls},C)$
$\displaystyle=-\log\bigg{(}\frac{\exp(h_{cls}*c_{i})}{\exp(h_{cls}*c_{i})+\sum\exp(h_{cls}*c_{i}^{\prime})}\bigg{)}$
(2)

Finally, we define the total loss as $L=\alpha L_{mlm}+\beta L_{nur}$. It’s
worth noting that we use the NUR task for effective vision and dialog fusion
rather than NSP or ITM tasks. Because these tasks are not only weaker than
seq2seq or bidirectional LM but also computationally expensive [19].

##  3 EXPERIMENT

Datasets. We train VU-BERT on the large visual dialog datasets VisDial v0.9
and v1.0. Specifically, the v1.0 contains 123,287 images for training, 2,064
images for validation and 8,000 images for test, which is a re-organization of
train and validation splits from v0.9 to form the new train. Each image for
training and validation is associated with one caption and 10 question-answer
pairs. For each question, it is paired with a list of 100 answer candidates,
one of which is regarded as the correct answer. In addition, visdial v1.0
provides dense annotations that specify real-valued relevance scores for the
100 answer candidates.

Baseline. We compare VU-BERT against a variety of single and ensemble
baselines in VD tasks, including FGA [9], VD-BERT [11], LF [12], HRE [12],
DVAN [13], DualVD [14], VisDial-BERT [15] and CAG [21], which yield prominent
performance from the leadboard or publication.

Experimental Setup. We use $BERT_{BASE}$ as the backbone, which consists of 12
Transformer blocks, each with 12 attention heads and a hidden state dimension
of 768. We set dropout probability to 0.1 and the learning rate with Adam
optimizer to 3e-4 and use GELU as activation function. In our experiments, we
split each $224\times 224$ image into a $7\times 7$ grid of image patches,
where each patch is $32\times 32$. We first train our model for 20 epochs with
a batchsize of 32 on a cluster of 8 Tesla V100 with 16G memory using MLM and
NUR losses.

We evaluate VU-BERT on VisDial v0.9 and v1.0 by using the ranking metrics
including Recall@K ($K\in\\{1,5,10\\})$: the percentage of the ground truth
answer option in top-k ranked responses, Mean: the average rank of the ground
truth answer option, and MRR (Mean Reciprocal Rank): reciprocal rank of the
ground truth answer option. Due to the acquisition of dense annotation, we add
the NDCG (Normalized Discounted Cumulative Gain), which penalizes the lower
rank of answers with high relevance, to evaluate each answer candidate for
v1.0.

##  4 Results and Analyses

Table 2: Performance of discriminative models on VisDial v1.0 dataset. Higher is better for MRR, recall@k and NDCG, while lower is better for Mean. Model | MRR$\uparrow$ | R@1$\uparrow$ | R@5$\uparrow$ | R@10$\uparrow$ | Mean$\downarrow$ | NDCG$\uparrow$
---|---|---|---|---|---|---
NMN | 0.5880 | 44.15 | 76.88 | 86.88 | 4.81 | 0.5810
GNN | 0.6137 | 47.33 | 77.98 | 87.83 | 4.57 | 0.5282
DualVD | 0.6323 | 49.25 | 80.23 | 89.70 | 4.11 | 0.5632
CAG | 0.6349 | 49.85 | 80.63 | 90.15 | 4.11 | 0.5664
LF | 0.5542 | 40.95 | 72.45 | 82.83 | 5.95 | 0.4531
FGA | 0.6370 | 49.58 | 80.97 | 88.55 | 4.51 | 0.5210
VD-BERT | 0.5117 | 38.90 | 62.82 | 77.98 | 6.69 | 0.7454
VisDial-BERT | 0.5074 | 37.95 | 64.13 | 80.00 | 6.53 | 0.7447
VU-BERT | 0.4909 | 33.60 | 67.20 | 81.60 | 6.12 | 0.7287
Table 3: The parameters for image extraction models Image model | VGG | Fast R-CNN | Patch
---|---|---|---
Method |  | LF;HRE;
---
NMN; GNN;
DVAN; FGA
| DualVD; CAG;
---
VD-BERT;
VisDial-BERT
VU-BERT
# Parameters | 138.36M | 41.76M | 0.59M

###  4.1 Results

On VisDial v0.9. We show both discriminative and generative results on v0.9
with the state-of-the-art baselines from previous published results and
leaderboard in Table 1. VU-BERT yields comparable results with the state of
the art in the discriminative and generative settings, for which the MRR is
0.6333 and 0.5403 respectively. This validates the effectiveness of our model
using a unified Transformer encoder.

on VisDial v1.0. As shown in Tabel 2, we compare our results on VisDial v1.0
datasets. We can found that VU-BERT outperforms most of the model across
various metrics. Compared to VD-BERT, VU-BERT results in 4% improvement in
R@5, 3% improvement in R@10 and 4% reduction in Mean while in 2% reduction in
NDCG and in 2% reduction in MRR. This indicates that VU-BERT can achieve
comparable performance even without deep convolution.

Obviously, VU-BERT is a simpler unified model without any region features or
deep convolution operations. We commission the unified transformer module to
process image and text simultaneously in place of an additional image
extractor like Fast R-CNN, which leads to significant time and parameter
efficiency [4]. We notice that VU-BERT encodes candidate sets independently
rather than integrates each candidate response at the input layer. Therefore,
VU-BERT loses the early interaction between the input and the candidates,
which coincidentally agrees with the results of some metrics lower than VD-
BERT. In spite of a slight decline on metrics, VU-BERT can easily extend the
candidate sets to a large amount, which facilitates adaptation to other tasks.

Parameters. In order to present the lightweight VU-BERT, we compare the Patch
Projection with two image extraction models commonly used in baseline methods:
VGG and Fast R-CNN in Table 3. We can observe that the Patch Projection is
much lighter when capturing the image feature. Then we selected several
structurally similar models, which deal with text features by transformers, to
compare the number of parameters in Table 4. Obviously, VU-BERT is much
lighter than existing transformer-based models.

###  4.2 Ablation Study

In Fig 2, we represent how dialog history affects VD tasks. VU-BERT is trained
with different dialog turns with image input and test in v1.0 validation set
according to the dense annotation. We found that the metrics with only one
turn or no dialog history is lower than other multi-turn dialog history. That
is to say directly employing the models specific to the VQA and VC for VD
tasks will not achieve better results. Also, more does not necessarily mean
better for the dialog history. When $\text{turns}>7$, the metric values of
Recall@K, MRR and NDCG gradually become flat, which may be caused by the
limitation of input length. Besides, we compare the performance with the image
input or not. As shown in Fig 2, the result of $\text{vis}=0$ is obviously
lower when the dialog turns are both 9, which proves that VU-BERT takes
advantage of multi-modal information to improve performance.

Table 4: The parameters of transformer-based architecture Model | Visdial-Bert | VD-Bert | VU-Bert
---|---|---|---
# Parameters | 315.26M | 174.51M | 133.94M
Fig. 2: Ablation study. $\text{vis}=0$ denotes the input includes 9 dialog
turns with no image while others denote the different turns with images.

##  5 conclusion

In this paper, we present a simple visual dialog architecture VU-BERT with the
effective fusion of two modalities. VU-BERT is comparable to competitors on VD
tasks which are heavily equipped with convolution networks like Faster R-CNN.
Besides, we explore the effect of dialog turns and found VQA and VC models are
not suitable applied directly in VD tasks. This paper is an initial attempt to
explore patch embedding adapted to VD tasks and we ask for future work on VD
to focus more on the light model and deep interactions between images and
dialog history with no convolution.

##  6 acknowledgment

This paper is supported by the Key Research and Development Program of
Guangdong Province under grant No. 2021B0101400003 and the National Key
Research and Development Program of China under grant No. 2018YFB0204 403.
Corresponding author is Jianzong Wang from Ping An Technology (Shenzhen) Co.,
Ltd (jzwang@188.com).

## References

  * [1] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova,  “BERT: pre-training of deep bidirectional transformers for language understanding,”  in Proceedings of NAACL, 2019, pp. 4171–4186.
  * [2] Shijing Si, Rui Wang, Jedrek Wosik, Hao Zhang, David Dov, Guoyin Wang, and Lawrence Carin,  “Students need more attention: Bert-based attention model for small data with application to automatic patient message triage,”  in Machine Learning for Healthcare Conference. PMLR, 2020, pp. 436–456.
  * [3] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  “Deep residual learning for image recognition,”  in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [4] Wonjae Kim, Bokyung Son, and Ildoo Kim,  “Vilt: Vision-and-language transformer without convolution or region supervision,”  in International Conference on Machine Learning, 2021, pp. 5583–5594.
  * [5] Kelvin Xu, Jimmy Ba, Ryan Kiros, Kyunghyun Cho, Aaron Courville, Ruslan Salakhudinov, et al.,  “Show, attend and tell: Neural image caption generation with visual attention,”  in International conference on machine learning, 2015, pp. 2048–2057.
  * [6] Yongfei Liu, Bo Wan, Xiaodan Zhu, and Xuming He,  “Learning cross-modal context graph for visual grounding,”  in Proceedings of the AAAI, 2020, vol. 34, pp. 11645–11652.
  * [7] Akira Fukui, Dong Huk Park, Daylen Yang, Anna Rohrbach, Trevor Darrell, and Marcus Rohrbach,  “Multimodal compact bilinear pooling for visual question answering and visual grounding,”  in Proceedings of EMNLP, 2016, pp. 457–468.
  * [8] Zongwen Bai, Ying Li, Marcin Woźniak, Meili Zhou, and Di Li,  “Decomvqanet: Decomposing visual question answering deep network via tensor decomposition and regression,”  Pattern Recognition, vol. 110, pp. 107538, 2021.
  * [9] Idan Schwartz, Seunghak Yu, Tamir Hazan, and Alexander G Schwing,  “Factor graph attention,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2019, pp. 2039–2048.
  * [10] Yashar Deldjoo, Johanne R Trippas, and Hamed Zamani,  “Towards multi-modal conversational information seeking,”  in Proceedings of SIGIR, 2021, vol. 21.
  * [11] Yue Wang, Shafiq R. Joty, Michael R. Lyu, Irwin King, Caiming Xiong, and Steven C. H. Hoi,  “VD-BERT: A unified vision and dialog transformer with BERT,”  in Proceedings of EMNLP, 2020, pp. 3325–3338.
  * [12] Abhishek Das, Satwik Kottur, Khushi Gupta, Avi Singh, Deshraj Yadav, José MF Moura, et al.,  “Visual dialog,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2017, pp. 326–335.
  * [13] Dan Guo, Hui Wang, and Meng Wang,  “Dual visual attention network for visual dialog.,”  in IJCAI, 2019, pp. 4989–4995.
  * [14] Xiaoze Jiang, Jing Yu, Zengchang Qin, Yingying Zhuang, Xingxing Zhang, Yue Hu, and Qi Wu,  “Dualvd: An adaptive dual encoding model for deep visual understanding in visual dialogue,”  in Proceedings of the AAAI, 2020, vol. 34, pp. 11125–11132.
  * [15] Vishvak Murahari, Dhruv Batra, Devi Parikh, and Abhishek Das,  “Large-scale pretraining for visual dialog: A simple state-of-the-art baseline,”  in European Conference on Computer Vision. Springer, 2020, pp. 336–352.
  * [16] Shaoqing Ren, Kaiming He, Ross Girshick, and Jian Sun,  “Faster r-cnn: Towards real-time object detection with region proposal networks,”  Advances in neural information processing systems, vol. 28, pp. 91–99, 2015.
  * [17] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, et al.,  “An image is worth 16x16 words: Transformers for image recognition at scale,”  in 9th International Conference on Learning Representations, 2021\.
  * [18] Yonghui Wu, Mike Schuster, Zhifeng Chen, Quoc V Le, Mohammad Norouzi, Wolfgang Macherey, et al.,  “Google’s neural machine translation system: Bridging the gap between human and machine translation,”  arXiv preprint arXiv:1609.08144, 2016.
  * [19] Luowei Zhou, Hamid Palangi, Lei Zhang, Houdong Hu, Jason Corso, and Jianfeng Gao,  “Unified vision-language pre-training for image captioning and vqa,”  in Proceedings of the AAAI, 2020, vol. 34, pp. 13041–13049.
  * [20] Li Dong, Nan Yang, Wenhui Wang, Furu Wei, Xiaodong Liu, Yu Wang, et al.,  “Unified language model pre-training for natural language understanding and generation,”  in Advances in Neural Information Processing Systems, 2019, pp. 13042–13054.
  * [21] Dan Guo, Hui Wang, Hanwang Zhang, Zheng-Jun Zha, and Meng Wang,  “Iterative context-aware graph inference for visual dialog,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020, pp. 10055–10064.