Two Channel Filter Banks on  Arbitrary Graphs with Positive Semi Definite Variation Operators

Eduardo Pavez, Benjamin Girault, Antonio Ortega,  Philip A. Chou

E. Pavez and A. Ortega  are with the Department of Electrical and Computer Engineering,
University of Southern California, Los Angeles,
CA, USA. (e-mail<EMAIL_ADDRESS>ortega@sipi.usc.edu).
B. Girault is with Université de Rennes, ENSAI, CNRS, CREST-UMR 9194, Rennes, France. (e-mail: benjamin.girault@ensai.fr).
P. A. Chou is with Google Research, Seattle, Washington, USA. (e-mail:philchou@google.com).

We study the  design of  filter banks for signals defined on the nodes of graphs.
We propose novel two channel  filter banks, that  can be applied to  arbitrary graphs,  given a positive semi definite variation operator, while  using downsampling operators on arbitrary vertex partitions.
The  proposed filter banks also satisfy several desirable properties, including perfect reconstruction,  and  critical sampling, while having efficient implementations.  Our results generalize previous approaches only valid for the normalized Laplacian of bipartite graphs.
We consider  graph Fourier transforms (GFTs) given by the generalized eigenvectors of the variation operator. This  GFT basis is orthogonal in an alternative inner product space, which depends on the choices of downsampling sets and  variation operators.
We show that the spectral folding  property of the normalized Laplacian of bipartite graphs, at the core of bipartite filter bank theory, can be generalized for the proposed  GFT if the inner product matrix is chosen properly.
We  give a probabilistic interpretation to  the proposed  filter banks  using  Gaussian graphical models. We also study orthogonality properties of tree structured filter banks, and propose a vertex partition algorithms for downsampling. We show that the proposed filter banks can be implemented efficiently on  3D point clouds, with hundreds of thousands of points (nodes), while also improving the color signal representation quality over competing state of the art approaches.

two channel filter banks, graph fourier transform, graph signal, multiresolution representation

§ INTRODUCTION

Graphs are powerful tools to model unstructured data. Graph's nodes and edges represent objects of interest and  their similarity relations, respectively. A function on the nodes is called a graph signal, with examples including color attributes in images and  3D point clouds,  measurements obtained from   sensor networks,  or biomedical brain signals  [1]. Graph signal processing (GSP) focuses on developing new theories and algorithms for restoration, compression, and analysis of graph signals [2, 3, 4].

Filter banks and other multiresolution representations (MRR)    have been widely applied to signal processing  problems on regular grids [5, 6, 7, 8], which  has motivated their  generalization to graphs [2, 9].

Bipartite Filter Banks (BFB) [10, 11] are critically sampled two channel filter banks on bipartite graphs.
They are
constructed with low and high pass spectral
graph filters of the normalized graph Laplacian and  down/up sampling operators (see   <ref>).
BFBs satisfy several desirable properties such as:
(i)  perfect reconstruction (PR), (ii)  critical sampling (non redundant), (iii)  compact support (polynomial  filters), and  (iv)  energy preservation (orthogonality).
BFBs are closely related to  filter banks on regular grids, since both involve satisfying the same design conditions  (i)-(iv) and  specific design techniques  (e.g., Meyer and  Daubechies filters [5, 6])
can be adapted to provide filter designs in the graph domain [10, 11, 12, 13]. Given their efficient implementations, close relation with  filter banks on regular grids, and their
well-understood theoretical properties, BFBs   have found  numerous applications in  compression [14, 15, 16, 17], denoising [18] and signal analysis [19, 20, 21].
However, a major limitation of  BFBs  is that they are only valid for bipartite graphs, using either the normalized Laplacian or the adjacency matrix. In applications, the given graph is rarely bipartite, so that bipartite graph learning or bipartite approximation are required [22, 23, 24, 25, 26], which increases computational complexity significantly, and may hinder the  signal representation quality.
Another limitation of BFBs comes from their use of the normalized Laplacian
(with lowest frequency eigenvector entries  proportional to the node degree) instead of the combinatorial Laplacian (constant lowest frequency eigenvector
entries), which is usually a better choice for signal representations and
results in better energy compaction [27, 11, 28].

two channel filter bank implemented with spectral graph filters. Low and high pass analysis filters are denoted by $\Hm_0$ and $\Hm_1$ respectively.  $\Gm_0$ and $\Gm_1$ correspond to low and high pass synthesis filters.  $\Ac$ and $\Bc$ are the sampling sets of the low and high pass channels respectively.  $\xv$ and ${\yv}$ denote input and output signals respectively. $\av$ is the vector of approximation (low pass) coefficients, and $\dv$ is the vector of detail (high pass) coefficients.

For discrete time signals, many filter bank designs obey properties (i)–(iv) [5] (for example  JPEG2000 uses near orthogonal separable filter banks satisfying (i)–(iii) [29]).  These filter banks can be  designed in the frequency domain using the Discrete Fourier Transform (DFT), and  they  can be implemented efficiently by exploiting lifting structures, separable filtering and  the   Fast Fourier Transform (FFT).
In contrast,    fast Graph Fourier Transform (GFT) algorithms  are only available for very limited families of graphs  [30], while separable and lifting structures are only known for bipartite graphs [10, 11, 31],  thus preventing the use of these techniques for reduced complexity implementations on arbitrary graphs.
We are  motivated by 3D point cloud applications, where graphs can have  millions of nodes [32], thus preventing the use of high complexity filter banks, such as those that require eigendecomposition or  bipartite graph learning.  As a result, we seek   designs that can  leverage graph sparsity for implementation, for instance those using low degree polynomial graph filters.

In this paper we generalize the critically sampled filter bank design problem so that the bipartite graph and normalized Laplacian assumptions are no longer required to find a solution.

Instead, Generalized Filter Bank (GFB) solutions can be found for arbitrary graphs, for any positive semi-definite graph variation operator,  and using arbitrary
vertex partitions for downsampling,
while still satisfying properties (i)–(iv).  We show that GFBs have eigendecomposition free  implementations,  relying only on sparse matrix vector products, sparse linear solvers, and low complexity graph construction algorithms based on approximate $k$ nearest neighbors.

The key  innovation enabling GFBs is the adoption of a new inner product, which allows us to depart from the traditional Hilbert space (induced by the dot product) underlying  the majority of GSP methods[The use of non traditional Hilbert spaces has proven effective in various theoretical studies involving irregularly structured data, for example  graph  sparsification [33], semi-supervised and unsupervised learning [34],  compression of 3D point clouds [35, 36], and graph signal sampling [37].]. We build upon [28], in which  a graph is represented by a positive semi definite  variation operator $\Mm \succeq 0$   to measure signal smoothness, an  inner product   $\langle \xv, \yv \rangle_{\Qm} = \yv^{\top} \Qm \xv$, with $\Qm \succ 0$, and  the $(\Mm,\Qm)$ Graph Fourier Transform ($(\Mm,\Qm)$-GFT) given by the generalized eigenvectors of $\Mm$,  which form a $\Qm$-orthonormal basis.

While orthogonal solutions (i.e., leading to an $\Qm$-orthonormal basis for $\Qm = \Id$) are preferable in many scenarios (e.g., for coding applications), we will show that, for a given graph,  it is possible to find nearly-orthogonal solutions via appropriate choices of downsampling set and variation operator. Indeed, nearly-orthogonal solutions are often used in practice, if they provide other useful properties. As an example, orthogonal  filter banks cannot be constructed with finite impulse response linear phase filters (other than the Haar filters) [5] nearly-orthogonal bi-orthogonal filters have been used because of the advantages of symmetry for image coding applications  [29].

The BFB framework  is based on a spectral folding phenomenon  satisfied by the eigenvectors and eigenvalues of the normalized Laplacian  of bipartite graphs [38].
We prove a new spectral folding property for the $(\Mm,\Qm)$-GFT.

We show that for a given variation operator $\Mm$, and vertex partition for downsampling, there exists a unique inner product matrix $\Qm$ such that the $(\Mm,\Qm)$-GFT obeys the spectral folding property.  Based on this result,  we propose    perfect reconstruction, orthogonal and  biorthogonal filter banks using spectral graph filters  of the fundamental matrix $\Zm = \Qm^{-1}\Mm$. Our conditions in the graph frequency domain are exactly  those developed in  [10, 11] for the normalized Laplacian of bipartite graphs and therefore we can reuse any of previously proposed filter designs, including those in [10, 11]  as well as more recent improved designs, such as those in  [12, 13].
When the variation operator $\Mm$ is the normalized Laplacian, and the graph is bipartite, we recover the BFB framework.

In our preliminary version of this work [39], we first introduced the spectral folding property and two channel filter banks on arbitrary graphs (Sections <ref> and <ref>), while
an application of GFBs to  image compression has been recently published [27].
In this paper,
we further develop the GFB theory,
providing all proofs not given in [39], discussing the results in more depth and introducing several novel contributions, including:

* An analysis of tree structured iterated two channel GFBs, and corresponding frame bounds (<ref>).

* A probabilistic derivation of the  spectral folding property, and a probabilistic interpretation of   tree structured filter banks using Gaussian graphical models [This Gaussian  assumption is common in the GSP literature. It has been   used  for graph learning [40], graph sampling [41], while studying graph processes [42, 43], and compression  [44, 45, 46, 47, 48].] (<ref>).

* An algorithm to optimize the vertex partition for downsampling (<ref>).

§.§ Related  work

Ref.					 		    PR 		    O 		    Poly  	    CS		    AO 	    SO 	   GT 	    GO
[10] 		                    SGF 	    SGF 	    Bip 	    $\Lcb$

[11] 		                    SGF 	    SGF 	    Bip 	    $\Lcb$

[49, 50]                    SGF 	    SGF 	    $\Omega$ 	    $\Wm$

[13, 31]		                    SGF 	    SGF 	    Bip 	    $\Wm$
[51] 	                    SGF 	    SGF 	    Any 	    $\Lcb$

[52] 	                    SGF 	    pinv 	    Any 	    $\Mm$

[53]                    SGF     Interp     Any     Any

[54] 			                    SGF 	    Interp 	    Any 	    $\Mm$

[55]		                    SGF 	    SGF 	    Any 	    Any

Ours										            ($\Zm$)         SGF 	    SGF 	    Any 	    $\Mm$

Comparison of multiresolution representations for graph signals. Transform properties: perfect reconstruction (PR), orthogonal (O), polynomial filter (Poly), critical sampling (CS), analysis operator (AO), synthesis operator (SO), graph type (GT), and graph operator (GO). Analysis operators are implemented with spectral graph filters (SGF), while synthesis operators may be implemented via pseudo inverse (pinv), interpolation (Interp), or a SGF.

We discuss related filter banks and multi resolution representations on graphs from three perspectives:  graph topology, graph  matrix (variation operator), and  downsampling sets.

Graph topology.
Several existing   filter bank theories are applicable only to graphs with certain types of  topology, including  bipartite  [10, 11], $M$-block cyclic  [49, 50], circulant  [56], and  acyclic [57].
While  some of these graph structures appear naturally in some applications [58, 59],
in many cases graphs of interest do  not belong to any of these categories.
To apply BFBs to arbitrary graphs, one can decompose any graph as a sum of bipartite graphs and apply the filter bank in a separable manner [60]. Alternative approaches include bipartite approximation [22],  bipartite graph learning methods [23, 24, 25, 61, 26], graph oversampling [51], and vertex partition optimization [62].

Earlier filter bank designs for  arbitrary graphs  are difficult to invert (e.g., least squares reconstruction is needed) [63, 64, 65].
More recent multiresolution representations  fail to be simultaneously perfect reconstruction and orthogonal [52], while others, require full eigendecomposition [53, 66].

Existing approaches that are valid for arbitrary graphs have several disadvantages. On the one hand, approximation-based methods may reduce signal representation quality while also requiring additional computational resources (to select the best bipartite approximation) [22, 60, 23, 24].
On the other hand, methods that allow the original graph to be used can do so at
the expense of other desirable features, such as low complexity [53, 66], perfect reconstruction  or orthogonality [62, 52] (see <ref> for a comparison of some of these approaches).  Therefore,
a theoretical formulation leading to critically sampled filter banks for arbitrary graphs, without the aforementioned disadvantages,  can be an attractive alternative.

Graph operator.
The choice of graph operator,   such as the graph Laplacian or the adjacency matrix, whose non-zero pattern encodes the graph structure, is an important design decision  [67, 68]. Existing filter bank frameworks are built for specific types of graph operators because of their special algebraic or spectral properties [10, 49], but these choices may not  be suitable for a particular application. As an example, BFBs [10] use the normalized Laplacian, while the random walk Laplacian has been shown to achieve better coding performance [27]. More recent filter bank approaches [54, 55] can be applied to the larger class of positive semi-definite variation operators.  However,  they lack perfect reconstruction   [54] or  require computing a full eigendecomposition of the graph operator   [55], which significantly limits their  application to large graphs.
The proposed GFBs   are valid for positive semi definite variation operators, including graph Laplacians,  without compromising on other desirable properties.

For   discrete-time signals,  a downsampling by $2$ operator  keeps “every other sample” and discards the rest.  In graphs, there is no obvious notion of “every other vertex”, unless the graph is bipartite [60].  For graph filter banks,  vertex partitions of the node set can be chosen under various criteria [54, 62, 9].
Recently, [55, 69] used spectral domain sampling, and while this approach leads to an attractive theory, spectral sampling requires computing all eigenvectors and eigenvalues of the graph operator, which can have significant computation complexity.
[54] extended sampling theory of graph signals to filter banks, but  practical implementations of this framework cannot achieve perfect reconstruction. In this work,  we show that  any partition of the vertex set is a valid downsampling operator. We also propose strategies for optimally choosing these partitions in Section <ref>.

§.§ Paper outline

The rest of the paper is organized as follows.
In Sections <ref> and <ref> we introduce the fundamentals of GSP on arbitrary Hilbert spaces, and two channel filter banks on bipartite graphs, respectively. Our proposed inner product matrix $\Qm$, and GFBs are presented in   <ref>.   <ref> is dedicated to tree structured iterated filter banks.  In <ref> we give a probabilistic interpretation of  GFBs.
Section <ref> studies optimal vertex partitioning for downsampling.
We end this paper with  numerical results,    and conclusions in  Sections <ref>   and <ref> respectively. Proofs   can be found in the Appendix.

§ GSP IN ARBITRARY HILBERT SPACES

§.§ Notation

Scalars, vectors, and matrices are written in lower case regular, lower case bold, and upper case bold respectively (e.g.,  $a$, $\bv$, $\Cm$). Positive definite and semi-definite matrices are denoted by $\Am \succ 0$,  and $\Am \succeq 0$ respectively. We will denote by $\Cm_{\Ac \Bc}$, the sub-matrix of $\Cm$ whose rows and columns are indexed by the sets $\Ac$ and $\Bc$, respectively.  We  denote by  $\Vert \Am \Vert$, the spectral norm or largest singular value of a matrix $\Am$.

§.§ Graph signal processing

Consider a  graph $\Gc = (\Vc, \Ec)$ with vertex set $\Vc = \lbrace 1,2,\cdots,n \rbrace $, and edge set $\Ec \subset \Vc \times \Vc$.

A graph signal is a function $x: \Vc \rightarrow \mathbb{R}$, which can be represented by a vector $\xv = [x_1,\cdots,x_n]^{\top}$. The entry $x_i$ is the signal value at  vertex $i \in \Vc$.

The graph is equipped with a symmetric positive semi definite variation matrix $\Mm = (m_{ij})$, with sparsity pattern  determined by the edge set, that is, $m_{ij} = m_{ji}  = 0$ when $(i,j) \notin \Ec$ and $i \neq j$.

Throughout the paper we will assume that the graph is connected and thus  $\Mm$ is an irreducible matrix [70].
The graph signal variation is

\begin{equation}
\Delta(\xv) = \xv^{\top} \Mm \xv.
\end{equation}

Intuitively, given two signals of equal energy, the one with larger variation is the one with higher energy in the higher frequencies of the spectrum.
Typical choices of variation operator are the normalized,  combinatorial and generalized Laplacian matrices [71], which we define next.

The adjacency matrix is a non negative symmetric matrix $\Wm = (w_{ij})$, where $w_{ij} =0$, whenever $(i,j) \notin \Ec$. The degree of node $i$ is $d_i = \sum_{j \in \Vc }w_{ij}$, and the degree matrix is $\Dm = \diag(d_1,\cdots, d_n)$.
The combinatorial graph Laplacian (CGL) is $\Lm = \Dm - \Wm$, while the normalized graph Laplacian (NGL) is $\Lcb = \Dm^{-1/2}\Lm \Dm^{-1/2} = \Id - \Dm^{-1/2} \Wm \Dm^{-1/2}$.  A generalized graph Laplacian (GGL) is any positive semidefinite matrix with non positive off diagonal entries (see [71, 72, 40, 73] for examples).
Other variation operators are listed in  [28, 74].

Given a  positive definite matrix  $\Qm$, we consider the  following inner product between graph signals $\xv$ and $\yv$,

\begin{equation}
\langle \xv, \yv \rangle_{\Qm} = \yv^{\top} \Qm \xv,
\end{equation}

with   induced  $\Qm$-norm   $\Vert \xv \Vert_{\Qm} = \sqrt{\langle \xv, \xv \rangle_{\Qm}}$.

An orthonormal set in the $\Qm$ inner product, is a set $\lbrace \uv_k \rbrace_{k=1}^{n}$ that satisfies

\begin{equation}
\langle \uv_i, \uv_j \rangle_{\Qm} = \left\{ \begin{array}{cc}
1  & \mbox{if } i = j \\
0 & \mbox{if } i \neq j.
\end{array} \right.
\end{equation}

In matrix form this corresponds to $\Um^{\top}\Qm\Um = \Id$, where $\Um = [\uv_1,\cdots, \uv_n]$.
The $(\Mm,\Qm)$ Graph Fourier Transform ($(\Mm,\Qm)$-GFT)
is defined as the $\Qm$ orthonormal set that minimizes the graph signal variation [28]. The  $(\Mm,\Qm)$-GFT  can be obtained by solving the generalized eigendecomposition problem:

\begin{equation}
\Mm \uv_k = \lambda_k \Qm \uv_k.
\end{equation}

The matrix of generalized eigenvectors  is denoted by  $\Um = [\uv_1,\cdots, \uv_n]$, and the generalized eigenvalues can be arranged in the diagonal matrix $\Lambdam = \diag(\lambda_1,\cdots, \lambda_n)$, where   we also  assume that $\lambda_1 \leq \lambda_2,\cdots, \leq \lambda_n$. The set of generalized eigenvalues (or graph frequencies) is denoted by by $\sigma(\Mm,\Qm)$.

We can represent any graph signal $\xv$ in this basis as follows:

\begin{equation}\label{eq_MQGFT}
\xv = \sum_{i=1}^n \langle \xv, \uv_i \rangle_{\Qm} \uv_i = \Um \hat{\xv}.
\end{equation}

Here $\hat{\xv}$ is the $(\Mm,\Qm)$-GFT of $\xv$ [28].
In matrix form this  corresponds to $\hat{\xv} = \Um^{\top} \Qm \xv$.
The decomposition (<ref>) implements the inverse $(\Mm,\Qm)$-GFT, by reconstructing $\xv$ as a linear combination of the frequency components in $\hat{\xv}$.
This inverse transform is given by $\xv = \Um \hat{\xv}$, since $\Um \Um^{\top} \Qm = \Id$.

Consider a function $h: \sigma(\Mm,\Qm) \rightarrow\mathbb{R}$, then  a spectral graph filter   (SGF)  $\Hm$  is the matrix

\begin{equation}\label{eq_SGF_Definition}
\Hm = \Um h(\Lam) \Um^{\top} \Qm = h(\Zm),
\end{equation}

where  $h(\Lambdam) = \diag(h(\lambda_1), h(\lambda_2),\cdots, h(\lambda_n))$.  When $h(\lambda)=\lambda$, then
$\Hm=\Zm = \Qm^{-1} \Mm = \Um \Lam \Um^{\top} \Qm$ is called the fundamental matrix.

§.§ Vertex and frequency domain Hilbert spaces

Both a graph signal, $\xv$, and its $(\Mm,\Qm)$-GFT,  $\hat{\xv}$, can be represented by vectors in $\R^{n}$.
However, in our formulation they belong to different Hilbert spaces, namely,  $(\R^n,\langle \cdot, \cdot \rangle_\Qm)$ for $\xv$ and $(\R^n,\langle \cdot, \cdot \rangle_{\Id})$ for $\hat{\xv}$, where $\langle \cdot, \cdot \rangle_{\Id}$ is the traditional dot product.

Thus, we can view the $(\Mm,\Qm)$-GFT, $\Um^{\top}\Qm$, as a mapping from  $(\R^n,\langle \cdot, \cdot \rangle_\Qm)$ to $(\R^n,\langle \cdot, \cdot \rangle_{\Id})$, while the inverse $(\Mm,\Qm)$-GFT $\Um$ is a mapping from $(\R^n,\langle \cdot, \cdot \rangle_{\Id})$ to $(\R^n,\langle \cdot, \cdot \rangle_{\Qm})$.
Making explicit the Hilbert spaces involved in this mapping allows us to state the following property.

Let $\xv,\yv \in \R^{\Vc}$, and let $\hat{\xv}$ and $\hat{\yv}$ be their respective frequency domain representations, then

\begin{equation}\label{eq_ortho_GFT}
\langle \xv,\yv \rangle_{\Qm} = \langle \hat{\xv},\hat{\yv} \rangle_{\Id}.
\end{equation}

From    (<ref>) it is clear  that SGFs are functions from  $(\R^n,\langle \cdot, \cdot \rangle_{\Qm})$ to the same Hilbert space.
For the rest of the paper we will use the $\Qm$ inner product for graph signals, and the $\Id$ inner product for their  frequency domain representations.

§ CRITICALLY SAMPLED TWO CHANNEL FILTER BANKS

We formulate the two-channel filter bank design problem. In <ref>   we  define critically sampled two-channel filter banks without making assumptions about the graph, variation and downsampling operators.  Later in <ref>  we  review the solution obtained by  [10, 11] for bipartite graphs using the $(\Lcb,\Id)$-GFT.

§.§ Vertex domain conditions for arbitrary graphs

A two-channel filter bank is depicted in  <ref>, where $\Hm_0$ and $\Hm_1$ are the analysis filters, while  $\Gm_0$ and $\Gm_1$ are the synthesis filters.
Throughout the paper, we will consider critically sampled filter banks.

A two-channel filter bank is critically sampled if the downsampling sets $\Ac$ and $\Bc$ form a partition of $\Vc$.

Without loss of generality we assume that $\Ac = \lbrace 1, \cdots, \vert \Ac \vert \rbrace$, and $\Bc = \Vc \setminus \Ac$.  Downsampling a signal $\xv$ on a set $\Ac$ keeps the entries $\lbrace x_i: i \in \Ac \rbrace$, and discards the rest. A signal obtained by downsampling $\xv$ on $\Ac$ is denoted by

\begin{equation}
\xv_\Ac = \Sm_\Ac \xv,
\end{equation}

where $\Sm_\Ac = [\Id_\Ac, \mathbf{0}]$ is a $\vert \Ac \vert \times \vert \Vc \vert$ selection matrix.
The upsampling operator   $\Sm_\Ac^{\top}$ maps the signal back to $\Vc$ by filling the entries on $\Bc$ with zeroes. Downsampling followed by upsampling corresponds to

\begin{equation}
\begin{bmatrix}
\xv_\Ac^{\top} &
\mathbf{0}
\end{bmatrix}^{\top} = \Sm_\Ac^{\top} \xv_{\Ac} = \Sm_\Ac^{\top} \Sm_\Ac^{\phantom{\top}} \xv.
\end{equation}

The analysis operator $\Tm_a$ from <ref>  can be written

\begin{equation}\label{eq_Ta_analysis}
\Tm_a =  \Sm_\Ac^{\top} \Sm_\Ac^{\phantom{\top}} \Hm_0 + \Sm_\Bc^{\top} \Sm_\Bc^{\phantom{\top}} \Hm_1 = \begin{bmatrix}
\Sm_\Ac \Hm_0 \\ \Sm_\Bc \Hm_1
\end{bmatrix}.
\end{equation}

The output of the low pass and high pass channels are called approximation, $\av$, and detail, $\dv$, coefficients, respectively, and correspond to

\begin{equation}
\Tm_a \xv = \begin{bmatrix}
\av \\
\dv
\end{bmatrix}
= \begin{bmatrix}
\Sm_\Ac \Hm_0\xv \\ \Sm_\Bc \Hm_1 \xv
\end{bmatrix}.
\end{equation}

The synthesis operator takes $\Tm_a\xv$ as its input and has a similar expression

\begin{equation}\label{eq_synthesis_operator}
\Tm_s = \Gm_0 \Sm_\Ac^{\top} \Sm_\Ac^{\phantom{\top}}  + \Gm_1 \Sm_\Bc^{\top} \Sm_\Bc^{\phantom{\top}} = \begin{bmatrix}
\Gm_0 \Sm_\Ac^{\top} & \Gm_1 \Sm_\Bc^{\top}
\end{bmatrix}.
\end{equation}

Since the filterbank operates on graph signals, both the analysis and synthesis operators map signals from  $(\R^n,\langle \cdot, \cdot \rangle_\Qm)$ to $(\R^n,\langle \cdot, \cdot \rangle_\Qm)$ (see also <ref>).
Now that we have formally introduced two-channel  filter banks, we define their relevant properties  in terms of the operators $\Tm_a$ and $\Tm_s$.
A two-channel filter bank is perfect reconstruction (PR) if $\Tm_s \Tm_a \xv = \xv$, for all $\xv$.

Since we only consider critically sampled filter banks, the PR condition is reduced to $\Tm_s = \Tm_a^{-1}$, and a necessary and sufficient condition for perfect reconstruction  is given by

\begin{equation}\label{eq_PR_condition_sampling}
\Id = \Tm_s \Tm_a = \Gm_0 \Sm_\Ac^{\top} \Sm_\Ac \Hm_0  + \Gm_1 \Sm_\Bc^{\top} \Sm_\Bc  \Hm_1.  %\\
%&= \begin{bmatrix}
%  \Gm_0 \Sm_\Ac^{\top} & \Gm_1 \Sm_\Bc^{\top}
% \end{bmatrix}
%\Sm_\Ac \Hm_0 \\ \Sm_\Bc \Hm_1
\end{equation}

We also introduce  $\Qm$-orthogonal filter banks.

A  two channel filter bank is $\Qm$-orthogonal, if for every pair of graph signals  $\xv, \yv$,

\begin{equation}\label{eq_Q_ortho_FB}
\langle \xv, \yv \rangle_{\Qm} =  \langle \Tm_a \xv, \Tm_a \yv  \rangle_{\Qm}.
\end{equation}

In matrix form,  (<ref>)  is equivalent to $\Tm_a^{\top} \Qm \Tm_a = \Qm$.

Critically sampled $\Qm$-orthogonal filter banks are always perfect reconstruction. We can see this by noting that $\Tm_a^{-1}=\Tm_s = \Qm^{-1}\Tm_a^{\top} \Qm$. The synthesis operator $\Tm_s$ is also $\Qm$-orthogonal, since for any $\xv,\yv$

\begin{align}
\langle \Tm_s\xv, \Tm_s\yv \rangle_{\Qm} &= \langle \Tm_a \Tm_a^{-1}\Tm_s\xv,  \Tm_a \Tm_a^{-1}\Tm_s\yv \rangle_{\Qm} \\
&= \langle  \Tm_a^{-1}\Tm_s\xv, \Tm_a^{-1}\Tm_s\yv \rangle_{\Qm}  =\langle \xv, \yv \rangle_{\Qm}.
\end{align}

The definition of $\Qm$-orthogonality for filter banks is different than $\Qm$-orthogonality of the $(\Mm,\Qm)$-GFT. This is because the $(\Mm,\Qm)$-GFT maps graph signals from their vertex domain representation to their frequency domain representation, that is, $\Um^{\top}\Qm:~(\R^n,\langle \cdot, \cdot \rangle_\Qm) \rightarrow (\R^n,\langle \cdot, \cdot \rangle_\Id)$, which results in the  relationship of
<ref> (see [28]).

In contrast,  the analysis and synthesis operators of a PR critically sampled filter bank  stay in the same Hilbert space $(\R^n,\langle \cdot, \cdot \rangle_\Qm)$.

§.§ Spectral domain conditions for bipartite graphs

We seek two-channel filter banks that can provide a good representation for the signals of interest,
while having efficient   implementations on large (arbitrary) graphs. In practice, this can be achieved by designing SGFs $\Hm_i, \Gm_i$ that (i) have good frequency selectivity and
(ii) can be written as polynomials of the fundamental matrix $\Zm$.
In this section we review  the solutions from [10, 11] for bipartite graphs.  We start by defining bipartite graphs.

A graph $\Gc=(\Vc,\Ec)$  is  bipartite on  a partition of the node set $\Ac, \Bc$, if
for  all $(i,j)\in \Ec$, $i \in \Ac$ and $j \in \Bc$, or $i \in \Bc$ and $j \in \Ac$.

In bipartite graphs, only edges between  $\Ac$ and $\Bc$ are allowed. Thus,  the combinatorial Laplacian $\Lm$ and normalized Laplacian $\Lcb = \Dm^{-1/2}\Lm\Dm^{-1/2} = \Id - \tilde{\Wm} $
have the  form

\begin{equation}
\Lm = \begin{bmatrix}
\Dm_\Ac & -\Wm_{ \Ac\Bc} \\ -\Wm_{\Bc\Ac} & \Dm_\Bc
\end{bmatrix},~ \Lcb = \begin{bmatrix}
\Id_\Ac & -\tilde{\Wm}_{\Ac\Bc} \\ -\tilde{\Wm}_{\Bc\Ac} & \Id_\Bc
\end{bmatrix}.
\end{equation}

BFBs are implemented using spectral graph filters of  the $(\Lcb,\Id)$-GFT,
with $\Zm = \Lcb$, where for $i \in \{ 1,2 \}$,

\begin{equation}\label{eq_spec_graf_filt_normL}
\Hm_i = h_i(\Lcb),  \quad \Gm_i = g_i(\Lcb).
\end{equation}

The diagonal matrix $\Jm = \Sm^{\top}_{\Ac}\Sm_{\Ac}^{\phantom{\top}}  - \Sm^{\top}_{\Bc}\Sm_{\Bc}^{\phantom{\top}}$,
with entries

\begin{equation}\label{eq_f_setindicator}
\Jm_{i,i} = \left\{ \begin{array}{cc}
1  & \mbox{if } i \in \Ac \\
-1 & \mbox{if } i \in \Bc,
\end{array} \right.
\end{equation}

is used to establish the following property of the normalized Laplacian of bipartite graphs.

$\Lc \uv = \lambda \uv $ if and only if $\Jm \uv$ is also an eigenvector with eigenvalue $2 - \lambda$.

<ref> is the main property used by [10, 11]  to design perfect reconstruction and orthogonal filter banks.

For bipartite graphs with spectral graph filters given by (<ref>), a necessary and sufficient condition for PR is  that $ \forall \lambda \in \sigma(\Lcb,\Id)$

\begin{align}\label{eq_pr1}
g_0(\lambda)h_0(\lambda) +  g_1(\lambda)h_1(\lambda) &= 2, \\
h_1(\lambda)g_1(2-\lambda) - h_0(\lambda)g_0(2-\lambda) &= 0. \label{eq_pr2}
\end{align}

The proof  follows from using     $ \Sm_\Ac^{\top} \Sm_\Ac^{\phantom{\top}} = \frac{1}{2}(\Id +\Jm)$ and $\Sm_\Bc^{\top} \Sm_\Bc^{\phantom{\top}} = \frac{1}{2}(\Id -\Jm)$ in the PR condition  (<ref>), in combination with <ref> (see [10]).

<ref> can also be used to design  $\Id$-orthogonal  filter banks.

Under the same conditions of <ref>,  a  filter bank is $\Id$-orthogonal if and only if,
$ \forall \lambda \in \sigma(\Lcb,\Id)$

\begin{align}\label{eq_orthogonal_filters1}
h_0^2(\lambda) + h_1^2(\lambda) &= 2,\\
h_1(\lambda)h_1(2-\lambda) - h_0(\lambda)h_0(2-\lambda) &= 0.\label{eq_orthogonal_filters2}
\end{align}

An advantage of $\Id$-orthogonal filter banks, is that their design can be significantly simplified, by considering a single filter $h_0$, and setting $h_1(\lambda) = h_0(2 - \lambda)$, and $g_i = h_i$. This family  is known as Quadrature Mirror Filters (QMF) [5], and was extended to graphs in [10].

§ MAIN RESULTS

The proofs of the PR and orthogonal conditions for BFBs rely on the spectral folding property (<ref>), which only holds for the normalized Laplacian of bipartite graphs.
In this section we show that <ref> can be generalized to other graphs and variation operators if  we consider a different Hilbert space equipped with a particular type of   $\Qm$ inner product. Later, we apply this result to obtain Generalized Filter Banks on arbitrary graphs.

§.§ Lazy filter bank on arbitrary graphs

Lazy two channel filter bank

Before presenting our main results, we start with an example that illustrates how to extend the “lazy” filter bank (<ref>) from   bipartite   to  arbitrary graphs. Consider the matrix

\begin{equation}\label{eq_Cmatrix}
\Cm = \begin{bmatrix}
\Id_{\Ac} & \Cm_{\Ac \Bc} \\
\Cm_{\Bc \Ac} & \Id_{\Bc}
\end{bmatrix},
\end{equation}

where $\Cm_{\Ac \Bc}$ and $\Cm_{\Bc \Ac}$ are arbitrary rectangular matrices of dimensions $\vert \Ac \vert \times \vert \Bc \vert$ and $\vert \Bc \vert \times \vert \Ac \vert$, respectively.  It is easy to verify that  the lazy filter bank is  PR for arbitrary $\Cm$, indeed

\begin{align}
\hat{\xv} &= (\Sm_{\Bc}^{\top}\Sm_{\Bc}\Cm  + (2\Id - \Cm)\Sm_{\Ac}^{\top}\Sm_{\Ac})\xv \\
\mathbf{0} & \mathbf{0} \\
\Cm_{\Bc \Ac} & \Id_{\Bc}
\end{bmatrix}
+ \begin{bmatrix}
\Id_{\Ac} & \mathbf{0} \\
-\Cm_{\Bc \Ac} & \mathbf{0}
\end{bmatrix} \right)\xv = \xv.
\end{align}

A special case of this result is when  $\Cm=\Lcb$ and the graph is bipartite. While the matrix $\Cm$ and the normalized Laplacian have a similar structure,  since the  submatrices $\Cm_{\Ac \Bc}$ and $\Cm_{\Bc \Ac}$ are not constrained to the edge weights of a bipartite graph, we can choose the matrix $\Cm$ as a function of the variation operator $\Mm$ of an arbitrary graph.  In the following subsections, we will show that the choice of $\Cm_{\Ac \Bc} = \Mm_{\Ac \Ac}^{-1}\Mm_{\Ac \Bc}$ and  $\Cm_{\Bc \Ac} = \Mm_{\Bc \Bc}^{-1}\Mm_{\Bc \Ac}$,  can be used to construct generalized filter banks on arbitrary graphs.

§.§ Spectral folding on arbitrary graphs

First we define spectral folding for generalized eigenvectors.

Given a graph $\Gc$ with variation operator $\Mm\succeq 0$, inner product $\Qm \succ 0$ and a partition $\Ac$, $\Bc = \Vc \setminus \Ac$. The $(\Mm,\Qm)$-GFT has the spectral folding property,
if for  any generalized eigenpair $(\uv,\lambda)$ ,  $(\Jm\uv,(2-\lambda))$ is also a generalized eigenpair, that is:

\begin{equation}
\Mm \uv = \lambda \Qm \uv \Leftrightarrow  \Mm \Jm\uv = (2-\lambda) \Qm \Jm\uv,
\end{equation}

where $\Jm = \Sm_\Ac^{\top} \Sm_\Ac^{\phantom{\top}}  - \Sm_\Bc^{\top} \Sm_\Bc^{\phantom{\top}}$, as in (<ref>).

The following theorem is the main result of this paper, and completely characterizes the spectral folding property.

Given a graph $\Gc$ with variation operator $\Mm\succeq 0$, inner product $\Qm$ and a partition $\Ac = \lbrace  1, \cdots, \vert \Ac \vert \rbrace$, $\Bc = \Vc \setminus \Ac$. The $(\Mm,\Qm)$-GFT has the spectral folding property if and only if  $\Qm$ is chosen as

\begin{equation}
\Qm=\begin{bmatrix}\label{eq_M_Q}
\Mm_{\Ac\Ac} & \mathbf{0} \\ \mathbf{0} & \Mm_{\Bc\Bc}
\end{bmatrix}.
\end{equation}

Note that  we can always permute the rows and columns of $\Mm$, so that $\Ac = \{ 1,2,\cdots,\vert \Ac \vert \}$, up to vertex re-indexing.

<ref>  generalizes <ref>,  demonstrating that the spectral folding property is not unique to the  $(\Lcb,\Id)$-GFT of bipartite graphs, and in fact,  it is always satisfied if the inner product matrix  $\Qm$ is chosen as in (<ref>).
Note also that while other recent work [28, 35, 36, 37] demonstrates the benefits of using alternative inner products for signal processing on irregular domains, our work is the first to show that non-diagonal choices of $\Qm$ are also useful.

Using (<ref>),  the fundamental matrix $\Zm=\Um \Lambdam \Um^{\top}\Qm$ is

\begin{equation}\label{eq_fundamental_matrix_folding}
\Zm = \Qm^{-1}\Mm = \begin{bmatrix}
\Id_{\Ac} & \Mm_{\Ac\Ac}^{-1} \Mm_{\Ac \Bc} \\
\Mm_{\Bc\Bc}^{-1} \Mm_{\Bc\Ac} &\Id_{\Bc}
\end{bmatrix},
\end{equation}

which has the same structure of (<ref>), and thus $\Zm$ is a valid choice to construct  “lazy” filter bank from <ref>. In the next sections we will show that because of the spectral folding property, $\Zm$ in (<ref>) allows us to design a more general family of  filter banks.

In Appendix <ref> we derive some properties of the generalized eigenvalues that are consistent
with well-known results about the normalized and combinatorial Laplacian matrices.
Later in  <ref> we show that  $\Qm$ in (<ref>), can be alternatively derived using  probabilistic arguments based on Gaussian graphical models.

§.§ Proof of Theorem <ref>

First we  prove that  (<ref>) implies the spectral folding property (<ref>).
Since $\Mm$ and $\Qm$ are positive semidefinite, and $\Qm$ is non singular, there is a full set of generalized eigenvectors. Let $\uv = [\uv_\Ac^{\top}, \uv_\Bc^{\top} ]^{\top}$, then  $\Mm \uv = \lambda \Qm \uv$ is equivalent to

\begin{align}
\Mm_{\Ac\Ac} \uv_\Ac + \Mm_{\Ac\Bc} \uv_\Bc &= \lambda \Qm_\Ac \uv_\Ac,\label{eq_block_genEV1}  \\
\Mm_{\Bc\Ac} \uv_\Ac + \Mm_{\Bc\Bc} \uv_\Bc &= \lambda \Qm_\Bc \uv_\Bc.\label{eq_block_genEV2}
\end{align}

Set $\mathbf{v} = \Jm\uv=[\uv_\Ac^{\top}, -\uv_\Bc^{\top} ]^{\top}$, and using  (<ref>) and (<ref>)

\begin{align}
\Mm \vv &=
\begin{bmatrix}
\Mm_{\Ac\Ac} \uv_\Ac  -\Mm_{\Ac\Bc} \uv_\Bc \\
\Mm_{\Bc\Ac} \uv_\Ac  -\Mm_{\Bc\Bc} \uv_\Bc
\end{bmatrix} \\
\begin{bmatrix}
2 \Mm_{\Ac\Ac} \uv_\Ac  -\lambda \Qm_{\Ac} \uv_\Ac \\
\lambda \Qm_\Bc \uv_\Bc  - 2\Mm_{\Bc\Bc} \uv_\Bc
\end{bmatrix}=(2-\lambda) \Qm \vv.
\end{align}

The second equality is implied by $\Mm_{\Ac\Ac} = \Qm_\Ac$, and $\Mm_{\Bc\Bc} = \Qm_\Bc$. To prove the other direction of the equivalence in <ref>, we  set $\gamma=  2-\lambda$, and repeat the same steps.

Now, we prove that spectral folding implies <ref>.
We first show that $\mathbf{Q}$ is  block diagonal.

\begin{equation}
\Qm =
\begin{bmatrix}
\Qm_{\Ac}    & \Qm_{\Ac\Bc} \\
\Qm_{\Bc\Ac} & \Qm_{\Bc}
\end{bmatrix}.
\end{equation}

Using the spectral folding assumption, and given two generalized eigenvectors $\mathbf{u}$ and $\mathbf{v}$ with unit $\Qm$ norm, we have that $\Jm\uv$ and $\Jm\vv$ are also generalized eigenvectors of unit $\Qm$ norm. In addition,   $\Jm\uv$ is orthogonal to $\Jm\vv$, and $\uv$ is orthogonal to $\vv$.
Therefore, given the full matrix of generalized eigenvectors $\Um$, we have:

\begin{equation}\label{eq_orthoplusorth}
\Um^\top\left(\Qm+\Jm\Qm\Jm\right)\Um=2\Id
\text{.}
\end{equation}

Replacing  the right hand side of (<ref>) with $2 \Um^{\top}\Qm \Um = 2\Id$,  after simplification  we obtain $\Jm\Qm\Jm=\Qm$, meaning that $\Qm_{\Ac\Bc}=\mathbf{0}$, and $\Qm_{\Bc\Ac}=\mathbf{0}$.

Next, we prove that $\Qm_\Ac=\Mm_{\Ac \Ac}$ and $\Qm_\Bc=\Mm_{\Bc \Bc}$.
Using the spectral folding property, and given $\uv$, a generalized eigenvector with eigenvalue $\lambda$, we obtain:

\begin{align}
\Mm(\Id+\Jm)\uv &= (\lambda\Qm+(2-\lambda)\Qm\Jm)\uv \label{eq_folding_MAUA}\\
\Mm(\Id-\Jm)\uv &= (\lambda\Qm-(2-\lambda)\Qm\Jm)\uv \label{eq_folding_MBUB}.
\end{align}

(<ref>) and (<ref>) imply
$2\Mm_{\Ac \Ac}\uv_\Ac = 2\Qm_\Ac\uv_\Ac$ and $2\Mm_{\Bc\Bc}\uv_\Bc = 2\Qm_\Bc\uv_\Bc$, respectively.
Gathering all these equations for all eigenvectors into matrix form, we obtain $\diag(\Mm_\Ac,\Mm_\Bc)\Um=\diag(\Qm_\Ac,\Qm_\Bc)\Um$.
Since $\Um$ is invertible, we obtain the desired property that $\Mm_{\Ac \Ac}=\Qm_\Ac$ and $\Mm_{\Bc \Bc}=\Qm_\Bc$.

§.§ Vertex domain filtering and the fundamental matrix

Traditionally,  SGFs have been implemented using polynomials of the variation operator (i.e., the fundamental matrix of the $(\Mm,\Id)$-GFT).

Polynomial SGFs have efficient implementations based on sparse matrix vector products, while also being  interpretable in terms of local vertex domain operations.

Since $\Qm$ is not diagonal, unless the graph is bipartite,   $\Zm$ given by (<ref>) can be much denser than $\Mm$, which raises the issues of interpretability and computational complexity of polynomial SGFs of $\Zm$.
In this section we will show that the matrix vector product $\Zm \xv$, the key routine used to implement polynomial graph filters,  can be described in terms of vertex domain operations involving  bipartite and disconnected graphs. Complexity of SGFs as a function of the vertex partition is studied in  <ref>.

For simplicity, we  consider the combinatorial Laplacian, so that  $\Mm = \Lm$. For a given vertex partition $\Ac$ and $\Bc$,  any graph can be decomposed as the sum of a bipartite graph and a disconnected graph with $2$ or more connected components. An example is depicted in <ref>.
These graphs obey the following identities

\begin{equation}
\Lm = \Lbip + \Lblock, \quad \Wm = \Wbip + \Wblock, \quad \Dm = \Dbip + \Dblock,
\end{equation}

where the super indices bi and bd refer to bipartite and block diagonal, respectively (see
<ref>). In addition, we have that

\begin{equation}
\Lblock = \Dblock - \Wblock,\quad \Lbip = \Dbip - \Wbip.
\end{equation}

Representation of a graph (left) as the sum of a bipartite graph (center) and a disconnected graph with $2$ or more connected components (right). Red and blue circles contain node sets $\Ac$ and $\Bc$, respectively.

Now we can state the following.

The fundamental matrix in (<ref>) is equal to

\begin{equation}
\Zm = \Id - \Pblock \Pbip,
\end{equation}

where   $\Pbip = (\Dbip)^{-1}\Wbip$, and $\Pblock = ((\Dbip)^{-1}(\Dbip + \Lblock))^{-1}=((\Dbip)^{-1} \Qm^{-1}$, are both right stochastic non negative matrices, with bipartite and block diagonal structure, respectively.

The proof   can be found in Appendix <ref>.
The matrix $\Pblock \Pbip$ is a two step smoothing operator, since $\yv=\Pbip \xv$ is given by

\begin{equation}
y_i = (\Pbip \xv)_i = \left\{ \begin{array}{cc}
\frac{1}{(\Dbip)_{ii}}\sum_{j \in \Bc}w_{ij}x_j & \mbox{if } i \in \Ac \\
\frac{1}{(\Dbip)_{ii}}\sum_{j \in \Ac}w_{ij}x_j & \mbox{if } i \in \Bc.
\end{array} \right.,
\end{equation}

where a smoothed signal is obtained through linear combinations of neighbors on the complement set.
Since $\Pblock$ is block diagonal and non negative, the second filtering step is also a low pass filter that   only uses connections within $\Ac$ or $\Bc$. Complexity of the operation $\Zm \xv$ is thus dominated by the product with $\Pblock$ (i.e., product by $\Qm^{-1}$). In <ref> we will show that the vertex partition can be chosen so that the product $\Qm^{-1}\xv$ can be implemented efficiently.

§.§ Generalized filter banks on arbitrary graphs

We use the spectral folding property (<ref>) and   <ref>  from the previous section to construct perfect reconstruction (<ref>)  and $\Qm$-orthogonal (<ref>)   filter banks. Our proofs (Appendices <ref> and  <ref>) follow closely the proofs of <ref> and <ref> for the bipartite case, with the main difference being the use of spectral graph filters based on the $(\Mm,\Qm)$-GFT.

Now we present our PR result.

Let $\Mm$ be a positive semidefinite variation operator and let $\Ac$, $\Bc$ be a vertex partition. Let $\Qm$ be the inner product matrix chosen following <ref>  so that the $(\Mm,\Qm)$-GFT has the spectral folding property.
A  two channel filter bank with SGFs

\begin{align}\label{eq_spec_graf_filt_Z}
\Hm_i &= h_i(\Zm) = \Um h_i(\Lam)\Um^{\top}\Qm, \textnormal{ and }\\
\Gm_i &= g_i(\Zm) = \Um g_i(\Lam)\Um^{\top}\Qm
\end{align}

is perfect reconstruction, if and only if  for all $\lambda \in \sigma(\Mm,\Qm)$,

\begin{align}\label{eq_pr1_Q}
g_0(\lambda)h_0(\lambda) +  g_1(\lambda)h_1(\lambda) &= 2, \\
h_1(\lambda)g_1(2-\lambda) - h_0(\lambda)g_0(2-\lambda) &= 0. \label{eq_pr2_Q}
\end{align}

Next we present our  result for  $\Qm$ orthogonal filter banks.

Under the same conditions of  <ref>,   a two channel filter bank is $\Qm$-orthogonal if and only if,
for all $\lambda \in \sigma(\Mm,\Qm)$,

\begin{align}\label{eq_orthogonal_filters1_Q}
h_0^2(\lambda) + h_1^2(\lambda) &= 2,\\
h_1(\lambda)h_1(2-\lambda) - h_0(\lambda)h_0(2-\lambda) &= 0.\label{eq_orthogonal_filters2_Q}
\end{align}

Clearly, the PR and $\Qm$-orthogonality conditions  are the same as those required for   BFBs, and therefore  existing filter designs, such as those from  [10, 11, 12, 13],  can be reused.

If  $h_0$, $h_1$, $g_0$, and $g_1$  are  polynomials,    their  corresponding SGFs are   polynomials of $\Zm$.
In the BFB design, since  $\Lcb$ is sparse and $\Zm = \Lcb$, if the filters are low degree polynomials,    they can be implemented using localized vertex domain operations via sparse matrix vector products.
For GFBs, polynomial filters are also desirable, since they can  be implemented on large graphs by exploiting graph sparsity (see <ref>).

It was shown in [11] that the  $\Id$-orthogonal filter bank design problem does not have non-trivial polynomial solutions, i.e., of degree lower than that of the minimal degree polynomial of the graph operator. Unfortunately, these arguments also apply to $\Qm$-orthogonal filter banks.
An alternative  is to approximate   $\Qm$-orthogonal filters  with low degree Chebyshev polynomials [65, 12], in which case the filter bank can   be implemented on large graphs, at the expense of  satisfying  PR and  $\Qm$-orthogonality conditions only approximately.

Another approach is to use non $\Qm$-orthogonal PR filters.

Biorthogonal filters  defined by

\begin{equation}\label{eq_bior}
h_0(\lambda) = g_1(2-\lambda), \quad
h_1(\lambda) = g_0(2-\lambda),
\end{equation}

are perfect reconstruction.

Biorthogonal filter banks    have low degree polynomial implementations, and can  be designed to be near orthogonal [11].

We can use  the calculations of   Appendix <ref>,  and following  the  reasoning from <cit.>, it is easy to see that the biorthogonal GFBs from [11]  are approximately $\Qm$-orthogonal, and satisfy

\begin{equation}
\alpha  \Vert \xv \Vert_{\Qm} \lesssim \Vert \Tm_{*}\xv \Vert_{\Qm} \lesssim \beta \Vert \xv \Vert_{\Qm},
\end{equation}

for any $\xv$, where $*$ can be $a$ or $s$,   $a \lesssim b$ means there is some $\epsilon \geq 0$ so that $a \leq b + \epsilon$,   and the constants $\alpha,\beta$ are given by $\alpha^2 =(1/2) \inf_{\lambda \in [0,2]} (h_0^{2}(\lambda) + h_1^{2}(\lambda))$, and $\beta^2 =(1/2) \sup_{\lambda \in [0,2]} (h_0^{2}(\lambda) + h_1^{2}(\lambda))$.

The zeroDC filter bank[11] was  proposed so that DC (constant) signals are mapped to the lowest graph frequency ($\lambda =0$). This is achieved by multiplying the input signal by $\Dm^{1/2}$ before applying the analysis filter bank, and multiplying by $\Dm^{-1/2}$ at the output of the synthesis filter bank. This ensures that a constant input signal has zero response in the high pass channel. [11] showed that biorthogonal zeroDC filter banks can be implemented with polynomials of the  random walk Laplacian of a bipartite graph. The zeroDC filter banks can be derived as a special case of our framework with the $(\Lm,\Dm)$-GFT, by noticing that for bipartite graphs with Laplacian $\Lm$, the inner product matrix from  <ref>  is given by $\Qm = \Dm$.

§ MULTI-CHANNEL FILTER BANK

In this section we study properties of  tree structured filter banks formed by  concatenating  two-channel filter banks.  For all resolution levels, the  graphs, variation operators and sampling sets are given and fixed.  An example with $L=3$ levels is depicted in   <ref>.

Tree structured analysis filter bank

§.§ Analysis and synthesis operators

We assume the input signal is at  resolution $L$, thus $\av_L = \xv$. The outputs of the low and high  pass channels at resolution $\ell<L$ are called approximation and detail coefficients, and are denoted by    $\av_{\ell}$, and $\dv_{\ell}$,    respectively.

The sampling sets obey $\Vc_{\ell} = \Ac_{\ell}$, and for $\ell<L$, $\Ac_{\ell+1} = \Ac_{\ell}\cup\Bc_{\ell}$. The graph at resolution $\ell$  is denoted by $\Gc_{\ell} = (\Vc_{\ell},\Ec_{\ell})$, and has variation operator $\Mm_{\ell}$ with corresponding inner product matrix $\Qm_{\ell}$, chosen so that the $(\Mm_{\ell}, \Qm_{\ell})$-GFT  has the spectral folding property.

The analysis equation at resolution $\ell$ is given by

\begin{equation}
\begin{bmatrix}
\av_{\ell}^{\top} &
\dv_{\ell}^{\top}
\end{bmatrix}^{\top} =
\Tm_{a,\ell} \av_{\ell+1}.
\end{equation}

The synthesis operator is  $\Tm_{s,\ell}$, and satisfies

\begin{equation}\label{eq_synthesis_iterated_1level}
\av_{\ell+1} = \Tm_{s,\ell}\begin{bmatrix}
\av_{\ell}^{\top} &
\dv_{\ell}^{\top}
\end{bmatrix}^{\top}=\Thetam_{0,\ell} \av_{\ell} + \Thetam_{1,\ell} \dv_{\ell},
\end{equation}

where we have defined the matrices

\begin{equation}
\Thetam_{0,\ell} =  \Gm_{0,\ell} \Sm^{\top}_{\Ac_{\ell}}, \quad \Thetam_{1,\ell}= \Gm_{1,\ell} \Sm^{\top}_{\Bc_{\ell}},
\end{equation}

with dimensions $\vert \Ac_{\ell+1} \vert \times \vert \Ac_{\ell} \vert$,  and $\vert \Ac_{\ell+1} \vert \times \vert \Bc_{\ell} \vert $,  respectively.

After applying (<ref>) recursively we get the synthesis equation

\begin{equation}
\xv=\Thetam_{0,L-1} \av_{L-1} + \Thetam_{1,L-1} \dv_{L-1}
% &=\Thetam_{0,L-1} \left ( \Thetam_{0,L-2} \av_{L-2} + \Thetam_{1,L-2} \dv_{L-2} \right)+ \Thetam_{1,L-1} \dv_{L-1}\\
= \Psim_{0}\av_0 + \sum_{k=0}^{L-1} \Phim_k\dv_k.
\end{equation}

The matrices $\Psim_{k}$ and $\Phim_{k}$ are given by

\begin{equation}
\Psim_k = \Thetam_{0,L-1} \Thetam_{0,L-2}\cdots \Thetam_{0,k}, \quad
\Phim_k = \Psim_{k+1} \Thetam_{1,k},
\end{equation}

for $k \in [L-1]$, and $\Psim_{L}=\Id$. By collecting the coefficients into a vector $\cv = [\av^{\top}_0,  \dv^{\top}_0,    \cdots,   \dv^{\top}_{L-1}]^{\top}$,  we can write

\begin{equation}\label{eq_synthesis_tree1}
\xv =
\begin{bmatrix}
\Psim_{0} & \Phim_{0} &  \cdots & \Phim_{L-1}
\end{bmatrix}
\cv= \Tcb_s \cv,
\end{equation}

where $\Tcb_s$ and $\Tcb_a = \Tcb_s^{-1}$ are the synthesis and analysis operators of the tree structured filter bank, respectively.

§.§ Frame bounds

Consider a compression application, where $\cv$ and $\hat{\cv}$ are the unquantized and quantized coefficients, respectively. We would like to write the reconstruction error

\begin{equation}\label{eq_compression_error}
\xv - \hat{\xv} = \Tcb_s(\cv - \hat{\cv})
\end{equation}

as a function of the quantization error $\cv - \hat{\cv}$. The challenge, is that the   operators $\Tcb_a$ and $\Tcb_s$ are not orthogonal (under any norm). This is true even if each two channel filter bank at resolution $\ell$ is $\Qm_{\ell}$ orthogonal, in fact

\begin{equation}\label{eq_parseval_1level}
\Vert \av_{\ell+1}  \Vert_{\Qm_{\ell+1}}^2 = \Vert \av_{\ell}  \Vert_{\Qm_{0,\ell}}^2 + \Vert \dv_{\ell}  \Vert_{\Qm_{1,\ell}}^2,
\end{equation}

where we have  defined

\begin{equation}\label{eq_sub_QEll}
\Qm_{0,\ell} = \Qm_{\ell+1}(\Ac_{\ell}, \Ac_{\ell}),\quad
\Qm_{1,\ell} = \Qm_{\ell+1}(\Bc_{\ell}, \Bc_{\ell}).
\end{equation}

Note that if we wanted to write $\Vert \av_{\ell}  \Vert_{\Qm_{0,\ell}}^2$ as the sum of the energies of $\av_{\ell-1}$ and $\dv_{\ell-1}$, in some norm, and apply recursively (<ref>),  we would need that

\begin{equation}\label{eq_condition_parseval_ver1}
\Qm_{0,\ell} = \Qm_{\ell},
\end{equation}

which is satisfied by the normalized Laplacians of bipartite graphs at all resolution levels (since $\Qm_{0,\ell} = \Qm_{\ell}=\Id$). However, (<ref>) does not

hold for other variation operators on arbitrary graphs. Fortunately, we can still establish a relation between the norms of the coefficients at different resolutions. We will present a frame bound that relates the norms of $\xv$ and its coefficients $\cv$, when the two channel filter bank at level $\ell$ is $\Qm_{\ell}$ orthogonal. The frame constants depend on the amount by which (<ref>) is violated.  First we recall a basic result that relates different  vector norms.

Let $\Rm_1$ and $\Rm_2$ be positive definite matrices, then

\begin{equation}\label{eq_equivalent_norms}
\alpha \Vert \xv \Vert^2_{\Rm_2} \leq  \Vert \xv \Vert^2_{\Rm_1} \leq \beta \Vert \xv \Vert^2_{\Rm_2},
\end{equation}

for any $\xv$, where

\begin{equation}
\alpha = \inf_{\xv \neq 0} \frac{\xv^{\top}\Rm_1 \xv}{\xv^{\top}\Rm_2 \xv},\quad  \beta = \sup_{\xv \neq 0} \frac{\xv^{\top}\Rm_1 \xv}{\xv^{\top}\Rm_2 \xv}.
\end{equation}

Equalities in (<ref>) are attained if and only if $\Rm_1 = c \Rm_2$ for a given constant $c>0$, in which case $\alpha = \beta=c$.
For each $\ell \geq 0$, we define

\begin{equation}\label{eq_alpha_ell_beta_ell}
\alpha_{\ell} = \inf_{\xv \neq 0} \frac{\xv^{\top} \Qm_{0,\ell}  \xv}{\xv^{\top} \Qm_{\ell}  \xv}, \quad \beta_{\ell}=\sup_{\xv \neq 0} \frac{\xv^{\top} \Qm_{0,\ell}  \xv}{\xv^{\top} \Qm_{\ell}  \xv}.
\end{equation}

Also,  we define, $A_{L-1}=B_{L-1}=1$, and for $0 \leq \ell < L-1$,

\begin{equation}
A_{\ell} = \prod_{k=\ell+1}^{L-1} \alpha_k, \quad  B_{\ell} = \prod_{k=\ell+1}^{L-1} \beta_k.
\end{equation}

We also consider the following inner product matrix

\begin{equation}\label{eq_Q_iterated_filterbank}
{\Qm} = \begin{bmatrix}
\Qm_{0,0} & \mathbf{0} & \cdots & \mathbf{0} \\
\mathbf{0} & \Qm_{1,0} & \ddots & \mathbf{0}\\
\vdots & \ddots &  \ddots &  \vdots \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \Qm_{1,L-1}
\end{bmatrix}.
\end{equation}

Now we present or main result of this section for the tree structured filter bank of <ref>. The proof is given in Appendix <ref>.
For a  $L$-level  tree structured filter bank,  composed of $\Qm_{\ell}$ orthogonal two channel filter banks, the synthesis operator $\Tcb_s$ satisfies

\begin{equation}\label{eq_frame}
C \Vert \cv \Vert_{\Qm}^2 \leq \Vert \Tcb_s \cv \Vert^2_{\Qm_L} \leq D  \Vert \cv\Vert_{\Qm}^2,
\end{equation}

for all $\xv = \Tcb_s \cv$,
with constants given by

\begin{equation}
C = \min_{0 \leq \ell \leq L-1} A_{\ell}, \quad D = \max_{0 \leq \ell \leq L-1} B_{\ell}.
\end{equation}

Now we can revisit   (<ref>) and bound the reconstruction error using <ref>, leading to $\Vert \xv - \hat{\xv} \Vert^2_{\Qm_L} \leq D\Vert \cv - \hat{\cv}\Vert^2_{\Qm}$, where $\Qm$ is given by  (<ref>).

We have established that the synthesis operator $\Tcb_s$ is a frame with constants $C$ and $D$ [5, 8].   A frame is called tight, if the constants in the upper and lower bound are equal. In this case, $C=D$ when (<ref>) is true for all $\ell \geq 0$. GFBs are tight framew when all the graphs $\Gc_{\ell}$ are bipartite, and  BFBs are used  for each $\ell \geq 0$, since we have that $A = B =1$, $\Qm_L = \Qm = \Id$, and $\Vert \Tcb_s \cv \Vert_{\Id} = \Vert \cv\Vert_{\Id}$.

The constants $C$, $D$, and the inner product $\Qm_L$ and $\Qm$ that appear in <ref> depend on the graphs, variation operators and downsampling sets chosen at each resolution level.   Depending on the application, various criteria can be used to optimize the graphs, graph signal variation and downsampling operators.  In [27] we explored the design of variation operator for an image compression application. In  <ref> we study the vertex partition design problem.

§ A PROBABILISTIC INTERPRETATION OF GFBS

Assume that $\xv$ is a zero mean Gaussian vector with covariance matrix $\Sigmam$, and select $\Mm =\Sigmam^{-1}$  as the variation operator.   We will show that
$\Zm$ given by (<ref>) is the matrix that produces high pass coefficients $\Zm \xv$ with the least norm, hence  $\Zm$ solves

\begin{equation}\label{eq_optimality_Z}
\min_{\Cm} \E[\Vert \Cm \xv \Vert^2_{\Id}], \textnormal{ s.t. } \Cm = \begin{bmatrix}
\Id_{\Ac} & \Cm_{\Ac \Bc} \\
\Cm_{\Bc \Ac} & \Id_{\Bc}
\end{bmatrix},
\end{equation}

where $\Ac$ and $\Bc = \Vc\setminus \Ac$ are fixed, and  $\Cm$ is the high pass filter of the lazy two channel  filter bank from <ref>.
Later we study properties of the tree structured iterated lazy filter bank, and provide an alternative probabilistic perspective on its orthogonality.

We can reduce (<ref>) to the equivalent problem

\begin{align}
&\min_{\Cm_{\Ac \Bc}, \Cm_{\Bc \Ac}} \E[\Vert \xv_{\Ac} + \Cm_{\Ac \Bc}\xv_{\Bc} \Vert_{\Id}^2] + \E[\Vert \xv_{\Bc} + \Cm_{\Bc \Ac}\xv_{\Ac} \Vert_{\Id}^2]  \\
&= \min_{\Cm_{\Ac \Bc}} \E[\Vert \xv_{\Ac} + \Cm_{\Ac \Bc}\xv_{\Bc} \Vert_{\Id}^2] + \min_{ \Cm_{\Bc \Ac}} \E[\Vert \xv_{\Bc} + \Cm_{\Bc \Ac}\xv_{\Ac} \Vert_{\Id}^2] \nonumber.
\end{align}

Minimization of  $\E[\Vert \xv_{\Ac} + \Cm_{\Ac \Bc}\xv_{\Bc} \Vert_{\Id}^2]$ corresponds to optimal linear prediction of $\xv_{\Ac}$ from $\xv_{\Bc}$, which fortunately, has a closed form solution for Gaussian distributions.

Given a vertex partition $\Ac, \Bc$,  if $\Mm = \Sigmam^{-1}$, then $\Zm=\Qm^{-1}\Mm$ from (<ref>) is the minimizer of (<ref>), and

\begin{equation}
\Zm \xv
= \xv - \begin{bmatrix}
\E[\xv_{\Ac} | \xv_{\Bc}] \\
\E[\xv_{\Bc} | \xv_{\Ac}]
\end{bmatrix}
= \xv - \begin{bmatrix}
\Sigmam_{\Ac \Bc} \Sigmam_{\Bc \Bc}^{-1}\xv_{\Bc} \\
\Sigmam_{\Bc \Ac} \Sigmam_{\Ac \Ac}^{-1}\xv_{\Ac}
\end{bmatrix}.
\end{equation}

The proof of this result is a direct consequence of

\begin{equation}
\Sigmam_{\Ac \Bc} \Sigmam_{\Bc \Bc}^{-1} = - \Mm_{\Ac \Ac}^{-1}    \Mm_{\Ac \Bc},
\Sigmam_{\Bc \Ac} \Sigmam_{\Ac \Ac}^{-1} = - \Mm_{\Bc \Bc}^{-1}    \Mm_{\Bc \Ac},
\end{equation}

which follows  from  $\Mm \Sigmam = \Id$.

From the  $(\Mm,\Qm)$-GFT perspective, $\Zm$ is a high pass filter (see <ref>).  <ref> illustrates that the filtering operation $\Zm \xv$ is the prediction error of an optimal  linear predictor  for the Gaussian distribution.

Now we take the lazy filter bank from <ref> and form a tree structured filter bank as in <ref>. Filters at resolution $\ell$ are given by $\Hm_{0,\ell} = \Id$, $\Hm_{1,\ell} = \Zm_{\ell}$, $\Gm_{0,\ell} = 2\Id - \Zm_\ell$, and $\Gm_{1,\ell} = \Id$. We will show that  <ref> implies that this tree structured  filter bank  produces sub-bands with uncorrelated coefficients.

Assume that the input signal  $\xv = \av_L$ is a zero mean Gaussian with covariance $\Sigmam = \Mm^{-1}$. Since the low pass channel corresponds to down-sampling $\av_{\ell+1}$ on the set $\Ac_\ell$, we have that for  $\ell<L$

\begin{equation}\label{eq_lazy_a_ell}
\av_{\ell} = \xv_{\Ac_\ell},
\end{equation}

is a $\vert \Ac_\ell \vert$-dimensional zero mean Gaussian vector with covariance matrix $\Sigmam_{\Ac_\ell, \Ac_\ell}$. Therefore, at resolution $\ell$ we will use the variation operator $\Mm_\ell =\Sigmam_{\Ac_\ell, \Ac_\ell}^{-1} $, which can be computed using Schur complements  [75]

\begin{equation}
\Mm_{\ell-1} = (\Mm_{\ell})_{\Ac_\ell \Ac_\ell} - (\Mm_{\ell})_{\Ac_\ell \Bc_\ell}((\Mm_{\ell})_{\Bc_\ell \Bc_\ell})^{-1} (\Mm_{\ell})_{\Bc_\ell \Ac_\ell}.
\end{equation}

The detail coefficients   are

\begin{equation}\label{eq_lazy_d_ell}
\dv_{\ell} = \xv_{\Bc_\ell} - {\Sigmam}_{\Bc_\ell \Ac_\ell}\Sigmam^{-1}_{\Ac_\ell \Ac_\ell} \xv_{\Ac_\ell}.
\end{equation}

The coefficient vector obtained after iterating the filter bank $L$ times is given by $\cv = [\av_0^{\top}, \dv_0^{\top}, \cdots, \dv_{L-1}^{\top}]^{\top} = \Tcb_a \xv$.

The inverse covariance matrix of $\cv$ is equal to

\begin{equation}
(\E[\cv \cv^{\top}])^{-1} = \begin{bmatrix}
\Mm_{0}  & \mathbf{0} & \cdots & \mathbf{0} \\
\mathbf{0} & \Qm_{1,0} & \ddots & \mathbf{0}\\
\vdots & \ddots &  \ddots &  \vdots \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \Qm_{1,L-1}
\end{bmatrix},
\end{equation}

where the matrices $\Qm_{1,\ell}$ are defined in (<ref>).

The proof  follows from <ref>, and can be found in Appendix <ref>.
The iterated lazy  filter bank produces uncorrelated coefficients at different resolutions.

§ VERTEX PARTITIONING

In this section we propose a vertex partitioning algorithm for two channel GFBs.
We  assume  $\Mm$ is given, and we wish to optimize  $\Ac$ and $\Bc\setminus \Vc$, as a function of $\Mm$.
We start by formulating this problem as a  optimization of the condition number of $\Qm$. We provide some theoretical justifications for this choice of objective based on complexity and stability of GFBs. We end the section with our proposed solution.

§.§ Problem formulation

Because of  <ref>,  we can optimize the inner product matrix  $\Qm$ as a function of the vertex partition. We propose

\begin{align}\label{eq_min_condition}
\min_{\fv: \fv_i^2=1} \kappa(\Qm)\quad\textnormal{ s.t. } &\Qm = \frac{1}{2}(\Mm + \diag(\fv) \Mm\diag(\fv))\\
&\vert \sum_{i \in \Vc} \fv_i \vert = \parity( \vert \Vc \vert) \nonumber,
\end{align}

where $\kappa(\Qm)= \Vert \Qm \Vert \Vert \Qm^{-1}\Vert$ is the condition number of $\Qm$ and $\Vert \cdot \Vert$ is the operator norm (largest singular value). The constraint $\fv_i^2 = 1$ ensures that the vector $\fv$ has $i$th entry equal to $1$ when $i \in \Ac$, and equal to $-1$ otherwise.
Note that   $ \Jm = \diag(\fv)$, where  $\Jm$ is defined in (<ref>), while  the constraint $\Qm = \frac{1}{2}(\Mm + \diag(\fv) \Mm\diag(\fv))$ ensures that the $(\Mm,\Qm)$-GFT has the spectral folding property.  The function $\parity(k)$ is equal to $0$ when $k$ is even, and $1$ otherwise, therefore the constraint $\vert \sum_{i \in \Vc} \fv_i \vert = \parity( \vert \Vc \vert )$ guarantees that the partition is balanced. When there is an even number of nodes,  $\sum_{i \in \Vc} \fv_i = 0$ and  $\vert \Ac \vert = \vert \Bc \vert$. When  the number of nodes is odd, $\vert \sum_{i \in \Vc} \fv_i \vert= 1$, and $\vert \vert \Ac \vert - \vert \Bc \vert \vert =1$. Balanced partitions are desired, so that the middle eigenvalue $\lambda=1$ does not correspond to a high dimensional subspace (see Appendix <ref> for more details).

Next we justify the use of  $\kappa(\Qm)$ in   (<ref>).

Small condition number improves complexity  of GFBs

For large graphs, performing eigendecomposition to implement SGFs is infeasible due to high computational complexity, thus   polynomial graph filters are used instead. For BFBs, a polynomial graph filter  of degree $d$ can be implemented by sparse matrix vector products (of the normalized Laplacian) with complexity $\Oc(d\vert \Ec \vert)$  [10, 12].

For arbitrary graphs,   polynomials of $\Zm$ require the computation of products of the form $\Zm \xv$. However, although $\Mm$ may be sparse, the fundamental matrix $\Zm = \Qm^{-1}\Mm$ is not.  Hence,  a naive implementation of the product $\Zm\xv$ has complexity $\Oc(n^2)$, while direct computation of $\Zm = \Qm^{-1}\Mm$ has complexity $\Oc(n^3)$ due to matrix inversion.

A more efficient approach is to decompose the product $\yv = \Zm \xv$   into: 1) computing a sparse matrix vector product $\wv = \Mm \xv$,
and 2) solving  a   system of equations $\Qm \yv = \wv$.
The second step is the computation bottleneck, since it requires solving a  linear system.

For large scale problems, positive semidefinite linear system are solved by iterative methods such as the Conjugate Gradient (CG), whose complexity per iteration is that of  a sparse matrix vector products with $\Qm$. The condition number of $\Qm$ is directly responsable for: 1)   the convergence rate of CG algorithms (and thus the total number of iterations), and 2) stability to perturbations of the system $\Qm \yv = \wv$.   From a computational perspective, the best case scenario occurs when  $\Qm$ is sparse and has a small condition number. In practice this is attained when $\Qm$ is  diagonal, and the graph is bipartite.

Small condition number implies stable signal reconstruction In <ref> we showed that tree structured filter banks obey a frame bound with constants $C$ and $D$. To make the frame bounds as tight as possible, we can try to design the vertex partition so that  (<ref>) is approximately satisfied. Using the definition of $\alpha_{\ell}$ and $\beta_{\ell}$ from (<ref>) we have that

\begin{equation}
1 \leq  {\beta_{\ell}}/{\alpha_{\ell}} \leq \kappa(\Qm_{\ell}) \kappa(\Qm_{0,\ell}).
\end{equation}

Note that $\kappa(\Qm_{0,\ell})$ depends on the higher resolution variation operator $\Mm_{\ell+1}$, and is independent of $\Ac_{\ell}$. Since $\kappa(\Qm_{\ell})$ depends on  $\Mm_{\ell}$ and the set $\Ac_{\ell}$, we  can make  $\beta_{\ell}/\alpha_{\ell}$  closer to $1$ by minimizing $\kappa(\Qm_{\ell})$.

§.§ Approximate solution for generalized Laplacians

Because of the non convex constraint  $\fv_i^2=1$, the problem (<ref>) is also  non convex. In this section we obtain an upper bound for $\kappa(\Qm)$, which allows us to pose an alternative, more tractable optimization problem.
For the rest of this section we  assume that $\Mm$ is a positive semi definite generalized graph Laplacian, that is,  $\Mm = \Pm - \Wm$, where $\Pm$ is a diagonal matrix and $\Wm$ is a non negative matrix with zero diagonal.
Since $\Qm$ is given by <ref>, we have that  $\Qm = \Pm - \Wblock \succ 0$,  where $\Wblock = (1/2)(\Wm + \diag(\fv)\Wm\diag(\fv))$  is the block diagonal part of $\Wm$. We can   bound  $\kappa(\Qm)$ with

\begin{equation}\label{eq_bound_kappa1}
\kappa(\Qm) \leq \kappa(\Pm) ({1+\rho(\Ac)})/({1-\rho(\Ac)}),
\end{equation}

where $\kappa(\Pm)$ is the condition number of $\Pm$, and $\rho(\Ac) = \Vert \Pm^{-1/2}\Wblock \Pm^{-1/2} \Vert$. A derivation of this bound is given in Appendix <ref>.
When the graph is bipartite,  $\Qm$ is diagonal, $\rho(\Ac) = 0$,  and  the bound is tight.

Since the right side of (<ref>) is decreasing with $\rho(\Ac)$, we can minimize  $\rho(\Ac)$ as a function of the vertex partition.

As a simplification,  we use the  sequence of bounds for the operator norm $\Vert\Am \Vert \leq \Vert \Am \Vert_F \leq \Vert \Am \Vert_1 = \sum_{i,j}\vert a_{ij} \vert$, which results in $\rho(\Ac) \leq \sum_{i,j}(1+\fv_i \fv_j) w_{ij}/\sqrt{p_i p_j}$. We propose solving instead:

\begin{equation}\label{eq_min_L1_norm}
\min_{\fv: \fv_i^2=1}  \sum_{i,j}(1+\fv_i \fv_j)\tilde{w}_{ij}   \textnormal{ s.t. }   \vert \sum_{i \in \Vc} \fv_i \vert = \parity( \vert \Vc \vert)
\end{equation}

where $\tilde{w}_{ij}$ is the $ij$ entry of  $\tilde{\Wm} = \Pm^{-1/2}\Wm\Pm^{-1/2} $.

Using  the identity $\sum_{i,j}\fv_i \fv_j\tilde{w}_{ij} = \fv^{\top} \tilde{\Wm} \fv = \sum_{i,j}\tilde{w}_{ij} - \fv^{\top}\tilde{\Lm}\fv$, we have that
(<ref>) is equivalent to

\begin{equation}\label{eq_maxcut_f}
\max_{\fv: \fv_i^2} \fv^{\top} \tilde{\Lm} \fv  \text{ s.t. }   \vert \sum_{i \in \Vc} \fv_i \vert = \parity( \vert \Vc \vert),
\end{equation}

where  $\tilde{\Dm} = \diag(\tilde{\Wm} \mathbf{1})$,  and $\tilde{\Lm} = \tilde{\Dm} - \tilde{\Wm}$.  (<ref>) is an instance of weighted maximum cut (WMC) [76], with an additional balanced partition constraint
It is well known that WMC is NP hard, thus we consider a  spectral partitioning approximation [9, 77], that computes $\tilde{\uv}_n = \argmax_{\Vert \uv \Vert=1} \uv^{\top}\tilde{\Lm}\uv$, and sets  $\fv = \sign(\tilde{\uv}_n - \tau \mathbf{1})$, and $\Ac = \lbrace i \in \Vc: \fv_i =1 \rbrace$.
The parameter $\tau \in \R$ can be tuned to ensure $\vert \sum_{i \in \Vc} \fv_i \vert = \parity( \vert \Vc \vert)$.
In the next section, we show that this algorithm can be efficiently implemented  for large graphs (since it only requires computing an eigenvector), while  improving the signal representation.

Decomposition of $8$ connected $4 \times 4$ grid graph (<ref>) as sum of bipartite (<ref>) and block diagonal graphs (<ref>). Fundamental matrices of $8$ connected $8 \times 8$ grid graph using $(\Lm,\Qm)$-GFT (<ref>) and $(\Lm,\Dm)$-GFT (<ref>) have similar vertex domain localization.

GFB Meyer L=3

SDS Meyer L=3

GFB CDF-9/7 L=3

SDS CDF-9/7 L=3

Low frequency atoms  on Bunny point cloud centered at node $2000$. (left) Proposed  generalized filter banks (GFB), and (right) filter banks based on spectral domain sampling (SDS).

Approximation of 3d point cloud attributes using GFB.

§ EXPERIMENTS

In this section we present  numerical results comparing the proposed GFBs with existing designs.  We show that: 1) GFB can be implemented efficiently on graphs with hundreds of thousands of nodes, 2) the basis functions of GFBs are localized in the vertex domain, and 3) GFBs with optimized vertex partitions

can be better than BFBs for practical applications, through a point cloud representation example.

§.§ Vertex domain localization

We start by visualizing the vertex domain behavior of the  fundamental matrix $\Zm$ and the basis functions of the GFB.  We show that although   $\Zm$ is (numerically)  dense, in practice $\Zm$ is approximately sparse and   polynomial filters of $\Zm$ are localized in the vertex domain. Using (<ref>) (see Appendix) and since $\rho(A) <1$, we have for any $m \geq 1$

\begin{equation}\label{eq_Q_inv_neuman}
\Qm^{-1} =  \Pm^{-\frac{1}{2}}\left(\sum_{k=0}^{m}(\Pm^{-\frac{1}{2}}\Wblock \Pm^{-\frac{1}{2}})^{k}\right)\Pm^{-\frac{1}{2}} +\Em.
\end{equation}

The error matrix has norm $\Vert \Em\Vert  =  \Oc(\rho(\Ac)^{m+1})$, so that, when  $\rho(\Ac)$ is small,  $\Qm^{-1}$ is well approximated by an $m$-degree polynomial of a sparse matrix.

Fundamental matrix:
A typical graph used in image processing applications is the $8$-connected grid depicted in <ref>. The vertex partition (red or blue) can be used to decompose the graph into (i) a bipartite $4$-connected grid (<ref>), containing all vertical and horizontal connections,  and (ii) a block diagonal graph  containing  the diagonal connections (see <ref>). We compare  the random walk Laplacian (<ref>), which has the same sparsity as the $8$-connected grid, and the fundamental matrix based on the red/blue partition, with $\Qm$ chosen according to <ref>.  The proposed fundamental matrix (<ref>) and  the random walk Laplacian (<ref>) have comparable sparsity levels. We also observe that  the non zero entries of $\Zm$ decay in magnitude with the (graph) distance between nodes,  in agreement with (<ref>).

Basis functions of GFB:
We implement the tree structured filter bank from <ref>, and display the basis functions, given by the columns of $\Tcb_s$ from (<ref>). We compare them  to those of  spectral domain sampling (SDS) filter banks [55], which are  close to the proposed GFB in terms of graph properties  (see <ref>). Since SDS filter banks require full eigendecomposition, we use the relatively small Bunny point cloud, which has $n = 2503$ nodes, in our comparison. For both filter banks we fix the  number of levels $L=3$ and consider two designs: bi-orthogonal CDF-9/7 filters, and orthogonal Meyer filters.
We use  nonnegative kernel regression (NNK) for graph construction [78],  initialized with k-nearest-neighbors (KNN) with parameter $K=20$  and inverse distance as edge weights. For GFB, we use the spectral partitioning algorithm from <ref>.

<ref> displays low frequency basis functions, centered at a single node. For visualization purposes, we normalize the  values of  basis functions, so  their maximum entry is $1$. For both filter designs, the basis functions of the GFB   are  more localized than their SDS   counterparts. The  GFB with CDF-9/7 filters is implemented using polynomials of $\Zm$, thus producing the most localized basis functions.

§.§ Representation of  3D point cloud attributes

3D point clouds  consist of list of point coordinates $\Vm = [\vv_i] \in \R^{n \times 3}$, and  attributes $\av \in \R^{n}$.  A graph is constructed where each point is assigned to a node and the edge weight between nodes $i$ and $j$ is $w_{ij} = 1/\Vert \vv_i - \vv_j\Vert$. A sparse edge set is obtained using KNN and NNK [78] graph construction algorithms.
We consider the first frame of the loot and longdress sequences of the 8iVFBv2 dataset [32], which have $n=784,142$ and $n=765,821$ points respectively.  We implement the tree structured GFB of <ref> with $L=7$ levels,   using bi-orthogonal   analysis filters $h_0(\lambda) = \frac{1}{2a_0}(2-\lambda)(1+\lambda)$, and   $h_1(\lambda) = a_0 \lambda$, with  gain $a_0 = 0.735$, and synthesis filters computed with (<ref>).
We compare GFB using random sampling (GFB RaSa), GFB using the optimal sampling strategy from <ref> (GFB OpSa), BFB with combinatorial Laplacian and random sampling (BFB L), and BFB with normalized Laplacian and random sampling (BFB nL).
We approximate color signals by  keeping the low pass coefficients $\av_0$  and the  high pass coefficients    $\dv_i$ for $0 \leq i \leq m-1$, with $m \in {1,\cdots, L}$, and zeroing out the rest.
For each filter bank, we tested various graph constructions and chose the one giving best performance. For GFB, the best results are obtained with NNK graphs, while for BFB the best results are obtained with KNN graphs.  <ref> shows the PSNR between the color signal and its approximation.

As expected, the BFB with the normalized Laplacian has the worst performance. This can be attributed to the non constant DC signal [11, 27]. By using the BFB with the $(\Lm,\Dm)$-GFT, which results in the random walk Laplacian, we recover the zeroDC filter bank from [11], which has a much improved performance. The proposed GFB with random sampling always outperforms the BFB, although by a small margin. When using optimal sampling with GFB, performance is again improved consistently across datasets.

§ CONCLUSION

This paper proposed two-channel filter banks   on arbitrary graphs, with positive semi definite variation operators, and arbitrary vertex partitions for downsampling. Previous designs were only valid for the normalized Laplacian of bipartite graphs,  because its eigenvalues and eigenvectors   obey a spectral folding property, that is essential to ensure perfect reconstruction. Our main contribution is showing that this property is satisfied by the generalized eigenvalues and eigenvectors of arbitrary graphs, if the inner product matrix is properly chosen. Based on this, we proposed generalized filter banks (GFB), implemented with spectral graph filters of the generalized eigenvectors. We also studied other theoretical aspects of  GFBs, including iterated tree structured filters, probabilistic interpretations, and optimal vertex partitioning for downsampling. Numerical results show that GFBs      basis functions are localized in the vertex domain, they can be efficinently  implemented on large graphs with hundreds of thousands of nodes, while   outperforming bipartite filter banks on a signal representation tasks.

Some open questions and directions for future work include: better understanding of spectral properties of  $(\Mm,\Qm)$-GFTs, specially for non diagonal $\Qm$,  improved filter designs for the fundamental matrix, and application dependent   vertex partition algorithms.

[1]
A. Ortega, P. Frossard, J. Kovačević, J. M. Moura, and
P. Vandergheynst, “Graph signal processing: Overview, challenges, and
applications,” Proceedings of the IEEE, vol. 106, no. 5, pp.
808–828, 2018.

[2]
D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The
emerging field of signal processing on graphs: Extending high-dimensional
data analysis to networks and other irregular domains,” IEEE signal
processing magazine, vol. 30, no. 3, pp. 83–98, 2013.

[3]
A. Sandryhaila and J. M. Moura, “Discrete signal processing on graphs,”
IEEE transactions on signal processing, vol. 61, no. 7, pp.
1644–1656, 2013.

[4]
A. Ortega, Introduction to Graph Signal Processing.1em plus 0.5em
minus 0.4emCambridge University Press, 2022.

[5]
M. Vetterli and J. Kovacevic, Wavelets and subband coding.1em
plus 0.5em minus 0.4emPrentice-hall, 1995.

[6]
I. Daubechies, Ten lectures on wavelets.1em plus 0.5em minus
0.4emSiam, 1992, vol. 61.

[7]
P. P. Vaidyanathan, Multirate systems and filter banks.1em plus
0.5em minus 0.4emPearson Education India, 2006.

[8]
S. Mallat, A wavelet tour of signal processing.1em plus 0.5em
minus 0.4emElsevier, 1999.

[9]
D. I. Shuman, “Localized spectral graph filter frames: A unifying framework,
survey of design considerations, and numerical comparison,” IEEE
Signal Processing Magazine, vol. 37, no. 6, pp. 43–63, 2020.

[10]
S. K. Narang and A. Ortega, “Perfect reconstruction two-channel wavelet filter
banks for graph structured data,” IEEE Transactions on Signal
Processing, vol. 60, no. 6, pp. 2786–2799, 2012.

[11]
——, “Compact support biorthogonal wavelet filterbanks for arbitrary
undirected graphs,” IEEE transactions on signal processing, vol. 61,
no. 19, pp. 4673–4685, 2013.

[12]
A. Sakiyama, K. Watanabe, and Y. Tanaka, “Spectral graph wavelets and filter
banks with low approximation error,” IEEE Transactions on Signal and
Information Processing over Networks, vol. 2, no. 3, pp. 230–245, 2016.

[13]
D. B. Tay and J. Zhang, “Techniques for constructing biorthogonal bipartite
graph filter banks,” IEEE Transactions on Signal Processing, vol. 63,
no. 21, pp. 5772–5783, 2015.

[14]
H. Q. Nguyen, P. A. Chou, and Y. Chen, “Compression of human body sequences
using graph wavelet filter banks,” in 2014 IEEE International
Conference on Acoustics, Speech and Signal Processing (ICASSP).1em
plus 0.5em minus 0.4emIEEE, 2014, pp. 6152–6156.

[15]
A. Anis, P. A. Chou, and A. Ortega, “Compression of dynamic 3d point clouds
using subdivisional meshes and graph wavelet transforms,” in 2016 IEEE
International Conference on Acoustics, Speech and Signal Processing
(ICASSP).1em plus 0.5em minus 0.4emIEEE, 2016, pp. 6360–6364.

[16]
J. Zeng, G. Cheung, Y.-H. Chao, I. Blanes, J. Serra-Sagristà, and
A. Ortega, “Hyperspectral image coding using graph wavelets,” in 2017
IEEE International Conference on Image Processing (ICIP).1em plus
0.5em minus 0.4emIEEE, 2017, pp. 1672–1676.

[17]
D. E. O. Tzamarias, K. Chow, I. Blanes, and J. Serra-Sagristà,
“Compression of hyperspectral scenes through integer-to-integer spectral
graph transforms,” Remote Sensing, vol. 11, no. 19, p. 2290, 2019.

[18]
Y. Iizuka and Y. Tanaka, “Depth map denoising using collaborative graph
wavelet shrinkage on connected image patches,” in 2014 IEEE
International Conference on Image Processing (ICIP).1em plus 0.5em
minus 0.4emIEEE, 2014, pp. 828–832.

[19]
M. Levorato, S. Narang, U. Mitra, and A. Ortega, “Reduced dimension policy
iteration for wireless network control via multiscale analysis,” in
2012 IEEE Global Communications Conference (GLOBECOM).1em plus
0.5em minus 0.4emIEEE, 2012, pp. 3886–3892.

[20]
V. K. Sharma, D. K. Srivastava, and P. Mathur, “Efficient image steganography
using graph signal processing,” IET Image Processing, vol. 12, no. 6,
pp. 1065–1071, 2018.

[21]
Y.-l. Qiao, Y. Zhao, and X.-y. Men, “Target recognition in sar images via
graph wavelet transform and 2dpca,” in Proceedings of the 2nd
International Conference on Image and Graphics Processing, 2019, pp. 3–7.

[22]
S. K. Narang and A. Ortega, “Local two-channel critically sampled filter-banks
on graphs,” in 2010 IEEE International Conference on Image
Processing.1em plus 0.5em minus 0.4emIEEE, 2010, pp. 333–336.

[23]
J. Zeng, G. Cheung, and A. Ortega, “Bipartite approximation for graph wavelet
signal decomposition,” IEEE Transactions on Signal Processing,
vol. 65, no. 20, pp. 5466–5480, 2017.

[24]
A. Jiang, J. Wan, Y. Tang, B. Ni, and Y. Zhu, “Admm-based bipartite graph
approximation,” in ICASSP 2019-2019 IEEE International Conference on
Acoustics, Speech and Signal Processing (ICASSP).1em plus 0.5em minus
0.4emIEEE, 2019, pp. 5421–5425.

[25]
E. Pavez, H. E. Egilmez, and A. Ortega, “Learning graphs with monotone
topology properties and multiple connected components,” IEEE
Transactions on Signal Processing, vol. 66, no. 9, pp. 2399–2413, 2018.

[26]
S. Kumar, J. Ying, J. d. M. Cardoso, and D. P. Palomar, “A unified framework
for structured graph learning via spectral constraints,” Journal of
Machine Learning Research, vol. 21, no. 22, pp. 1–60, 2020.

[27]
D. E. O. Tzamarias, E. Pavez, B. Girault, A. Ortega, I. Blanes, and
J. Serra-Sagristà, “Exploiting orthogonality and zero-dc tradeoffs in
biorthogonal graph filterbanks,” in 2021 IEEE International Conference
on Acoustics, Speech, and Signal Processing (ICASSP).1em plus 0.5em
minus 0.4emIEEE, 2021.

[28]
B. Girault, A. Ortega, and S. S. Narayanan, “Irregularity-aware graph fourier
transforms,” IEEE Transactions on Signal Processing, vol. 66, no. 21,
pp. 5746–5761, 2018.

[29]
D. S. Taubman and M. W. Marcellin, JPEG2000 - image compression
fundamentals, standards and practice.1em plus 0.5em minus 0.4emKluwer, 2002.

[30]
K.-S. Lu and A. Ortega, “Fast graph fourier transforms based on graph symmetry
and bipartition,” IEEE Transactions on Signal Processing, vol. 67,
no. 18, pp. 4855–4869, 2019.

[31]
D. B. Tay and A. Ortega, “Bipartite graph filter banks: Polyphase analysis and
generalization,” IEEE Transactions on Signal Processing, vol. 65,
no. 18, pp. 4833–4846, 2017.

[32]
E. d’Eon, B. Harrison, T. Myers, and P. A. Chou, “8i voxelized full bodies-a
voxelized point cloud dataset,” ISO/IEC JTC1/SC29 Joint WG11/WG1
(MPEG/JPEG) input document WG11M40059/WG1M74006, 2017.

[33]
D. A. Spielman and N. Srivastava, “Graph sparsification by effective
resistances,” SIAM Journal on Computing, vol. 40, no. 6, pp.
1913–1926, 2011.

[34]
M. Hein, J.-Y. Audibert, and U. v. Luxburg, “Graph laplacians and their
convergence on random neighborhood graphs,” Journal of Machine
Learning Research, vol. 8, no. Jun, pp. 1325–1368, 2007.

[35]
P. A. Chou, M. Koroteev, and M. Krivokuća, “A volumetric approach to point
cloud compression—part i: Attribute compression,” IEEE Transactions
on Image Processing, vol. 29, pp. 2203–2216, 2019.

[36]
M. Krivokuća, P. A. Chou, and M. Koroteev, “A volumetric approach to point
cloud compression–part ii: Geometry compression,” IEEE Transactions
on Image Processing, vol. 29, pp. 2217–2229, 2019.

[37]
B. Girault, A. Ortega, and S. S. Narayayan, “Graph vertex sampling with
arbitrary graph signal hilbert spaces,” in 2020 IEEE International
Conference on Acoustics, Speech and Signal Processing (ICASSP).1em
plus 0.5em minus 0.4emIEEE, 2020, pp. 5670–5674.

[38]
F. R. Chung and F. C. Graham, Spectral graph theory.1em plus
0.5em minus 0.4emAmerican Mathematical Soc., 1997, no. 92.

[39]
E. Pavez, B. Girault, A. Ortega, and P. A. Chou, “Spectral folding and
two-channel filter-banks on arbitrary graphs,” in 2021 IEEE
International Conference on Acoustics, Speech, and Signal Processing
(ICASSP).1em plus 0.5em minus 0.4emIEEE, 2021.

[40]
H. E. Egilmez, E. Pavez, and A. Ortega, “Graph learning from data under
laplacian and structural constraints,” IEEE Journal of Selected Topics
in Signal Processing, vol. 11, no. 6, pp. 825–841, 2017.

[41]
A. Gadde and A. Ortega, “A probabilistic interpretation of sampling theory of
graph signals,” in 2015 IEEE international conference on Acoustics,
Speech and Signal Processing (ICASSP).1em plus 0.5em minus
0.4emIEEE, 2015, pp. 3257–3261.

[42]
C. Zhang, D. Florêncio, and P. A. Chou, “Graph signal processing-a
probabilistic framework,” Microsoft Res., Redmond, WA, USA, Tech. Rep.
MSR-TR-2015-31, 2015.

[43]
A. G. Marques, S. Segarra, G. Leus, and A. Ribeiro, “Stationary graph
processes and spectral estimation,” IEEE Transactions on Signal
Processing, vol. 65, no. 22, pp. 5911–5926, 2017.

[44]
E. Pavez, H. E. Egilmez, Y. Wang, and A. Ortega, “Gtt: Graph template
transforms with applications to image coding,” in 2015 Picture Coding
Symposium (PCS).1em plus 0.5em minus 0.4emIEEE, 2015, pp.

[45]
M. Rizkallah, X. Su, T. Maugey, and C. Guillemot, “Graph-based transforms for
predictive light field compression based on super-pixels,” in 2018
IEEE International Conference on Acoustics, Speech and Signal Processing
(ICASSP).1em plus 0.5em minus 0.4emIEEE, 2018, pp. 1718–1722.

[46]
C. Zhang and D. Florêncio, “Analyzing the optimality of predictive
transform coding using graph-based models,” IEEE Signal Processing
Letters, vol. 20, no. 1, pp. 106–109, 2012.

[47]
H. E. Egilmez, Y.-H. Chao, and A. Ortega, “Graph-based transforms for video
coding,” IEEE Transactions on Image Processing, vol. 29, pp.
9330–9344, 2020.

[48]
R. L. de Queiroz and P. A. Chou, “Transform coding for point clouds using a
gaussian process model,” IEEE Transactions on Image Processing,
vol. 26, no. 7, pp. 3507–3517, 2017.

[49]
O. Teke and P. P. Vaidyanathan, “Extending classical multirate signal
processing theory to graphs—part i: Fundamentals,” IEEE Transactions
on Signal Processing, vol. 65, no. 2, pp. 409–422, 2016.

[50]
——, “Extending classical multirate signal processing theory to
graphs—part ii: M-channel filter banks,” IEEE Transactions on Signal
Processing, vol. 65, no. 2, pp. 423–437, 2016.

[51]
A. Sakiyama and Y. Tanaka, “Oversampled graph laplacian matrix for graph
filter banks,” IEEE Transactions on Signal Processing, vol. 62,
no. 24, pp. 6425–6437, 2014.

[52]
D. I. Shuman, M. J. Faraji, and P. Vandergheynst, “A multiscale pyramid
transform for graph signals,” IEEE Transactions on Signal Processing,
vol. 64, no. 8, pp. 2119–2134, 2015.

[53]
S. Chen, R. Varma, A. Sandryhaila, and J. Kovačević, “Discrete
signal processing on graphs: Sampling theory,” IEEE transactions on
signal processing, vol. 63, no. 24, pp. 6510–6523, 2015.

[54]
S. Li, Y. Jin, and D. I. Shuman, “Scalable $ m $-channel critically sampled
filter banks for graph signals,” IEEE Transactions on Signal
Processing, vol. 67, no. 15, pp. 3954–3969, 2019.

[55]
A. Sakiyama, K. Watanabe, Y. Tanaka, and A. Ortega, “Two-channel critically
sampled graph filter banks with spectral domain sampling,” IEEE
Transactions on Signal Processing, vol. 67, no. 6, pp. 1447–1460, 2019.

[56]
M. S. Kotzagiannidis and P. L. Dragotti, “Splines and wavelets on circulant
graphs,” Applied and Computational Harmonic Analysis, vol. 47, no. 2,
pp. 481–515, 2019.

[57]
V. N. Ekambaram, G. C. Fanti, B. Ayazifar, and K. Ramchandran, “Spline-like
wavelet filterbanks for multiresolution analysis of graph-structured data,”
IEEE Transactions on Signal and Information Processing over Networks,
vol. 1, no. 4, pp. 268–278, 2015.

[58]
G. Cheung, E. Magli, Y. Tanaka, and M. K. Ng, “Graph spectral image
processing,” Proceedings of the IEEE, vol. 106, no. 5, pp. 907–930,

[59]
J.-Y. Kao, A. Ortega, and S. S. Narayanan, “Graph-based approach for motion
capture data representation and analysis,” in 2014 IEEE International
Conference on Image Processing (ICIP).1em plus 0.5em minus
0.4emIEEE, 2014, pp. 2061–2065.

[60]
S. K. Narang, Y. H. Chao, and A. Ortega, “Graph-wavelet filterbanks for
edge-aware image processing,” in 2012 IEEE Statistical Signal
Processing Workshop (SSP).1em plus 0.5em minus 0.4emIEEE,
2012, pp. 141–144.

[61]
K.-S. Lu, E. Pavez, and A. Ortega, “On learning laplacians of tree structured
graphs,” in 2018 IEEE Data Science Workshop (DSW).1em plus
0.5em minus 0.4emIEEE, 2018, pp. 205–209.

[62]
A. Anis and A. Ortega, “Critical sampling for wavelet filterbanks on arbitrary
graphs,” in 2017 IEEE International Conference on Acoustics, Speech
and Signal Processing (ICASSP).1em plus 0.5em minus 0.4emIEEE,
2017, pp. 3889–3893.

[63]
M. Crovella and E. Kolaczyk, “Graph wavelets for spatial traffic analysis,”
in IEEE INFOCOM 2003. Twenty-second Annual Joint Conference of the IEEE
Computer and Communications Societies (IEEE Cat. No. 03CH37428),
vol. 3.1em plus 0.5em minus 0.4emIEEE, 2003, pp. 1848–1857.

[64]
R. R. Coifman and M. Maggioni, “Diffusion wavelets,” Applied and
Computational Harmonic Analysis, vol. 21, no. 1, pp. 53–94, 2006.

[65]
D. K. Hammond, P. Vandergheynst, and R. Gribonval, “Wavelets on graphs via
spectral graph theory,” Applied and Computational Harmonic Analysis,
vol. 30, no. 2, pp. 129–150, 2011.

[66]
A. Cloninger, H. Li, and N. Saito, “Natural graph wavelet packet
dictionaries,” Journal of Fourier Analysis and Applications, vol. 27,
no. 3, pp. 1–33, 2021.

[67]
X. Dong, D. Thanou, M. Rabbat, and P. Frossard, “Learning graphs from data: A
signal representation perspective,” IEEE Signal Processing Magazine,
vol. 36, no. 3, pp. 44–63, 2019.

[68]
G. Mateos, S. Segarra, A. G. Marques, and A. Ribeiro, “Connecting the dots:
Identifying network structure via graph signal processing,” IEEE
Signal Processing Magazine, vol. 36, no. 3, pp. 16–43, 2019.

[69]
A. Miraki, H. Saeedi-Sourck, N. Marchetti, and A. Farhang, “Spectral domain
spline graph filter bank,” IEEE Signal Processing Letters, vol. 28,
pp. 469–473, 2021.

[70]
R. A. Horn and C. R. Johnson, Matrix analysis.1em plus 0.5em
minus 0.4emCambridge university press, 2012.

[71]
T. Biyikoglu, J. Leydold, and P. F. Stadler, Laplacian eigenvectors of
graphs: Perron-Frobenius and Faber-Krahn type theorems.1em plus 0.5em
minus 0.4emSpringer, 2007.

[72]
S. Kurras, U. Luxburg, and G. Blanchard, “The f-Adjusted Graph Laplacian: a
Diagonal Modification with a Geometric Interpretation,” in
Proceedings of the 31st International Conference on Machine
Learning, vol. 32, no. 2, Bejing, China, 22–24 Jun 2014, pp. 1530–1538.

[73]
E. Pavez, B. Girault, A. Ortega, and P. A. Chou, “Region adaptive graph
Fourier transform for 3D point clouds,” in 2020 IEEE International
Conference on Image Processing (ICIP).1em plus 0.5em minus
0.4emIEEE, 2020.

[74]
A. Anis, A. Gadde, and A. Ortega, “Efficient sampling set selection for
bandlimited graph signals using graph spectral proxies,” IEEE
Transactions on Signal Processing, vol. 64, no. 14, pp. 3775–3789, 2016.

[75]
F. Dorfler and F. Bullo, “Kron reduction of graphs with applications to
electrical networks,” IEEE Transactions on Circuits and Systems I:
Regular Papers, vol. 60, no. 1, pp. 150–163, 2012.

[76]
M. X. Goemans and D. P. Williamson, “Improved approximation algorithms for
maximum cut and satisfiability problems using semidefinite programming,”
Journal of the ACM (JACM), vol. 42, no. 6, pp. 1115–1145, 1995.

[77]
B. Aspvall and J. R. Gilbert, “Graph coloring using eigenvalue
decomposition,” SIAM Journal on Algebraic Discrete Methods, vol. 5,
no. 4, pp. 526–538, 1984.

[78]
S. Shekkizhar and A. Ortega, “Graph construction from data by non-negative
kernel regression,” in ICASSP 2020-2020 IEEE International Conference
on Acoustics, Speech and Signal Processing (ICASSP).1em plus 0.5em
minus 0.4emIEEE, 2020, pp. 3892–3896.

§.§ spectral folding property for adjacency

If a matrix has a bipartite sparsity pattern, then it has symmetric spectrum,
that is if $(\uv,\lambda)$ is an eigenvector-eigenvalue, then $\vv = [\uv_A^{\top} , -\uv_B^{\top}]^{\top}$ is also en eigenvector with eigenvalue $-\lambda$.
The above property also holds for generalized eigenvectors in any inner product.

§.§ Idea 1

consider the adjacency

\begin{equation}
\hat{\Bm}=  \Qm^{-1/2}\Bm \Qm^{-1/2}.
\end{equation}

This is like the normalized adjacency of the normalized laplacian. $\hat{\Bm}$ has the properties

* Symmetric, so orthogonal eigenvectors

* Possibly dense.

* $\Qm^{1/2}\mathbf{1}$ is an eigenvector with eigenvalue $1$.

* Eigenvalues in $[-1,1]$

§.§ Idea 2

We can look at the matrix

\begin{equation}
\Id - \Zm = \Qm^{-1}\Bm  = \tilde{\Bm}
\end{equation}

as the new adjacency matrix. $\tilde{\Bm}$ is right stochastic, so we can look at it as an alternative  to the degree normalization  when the graph is non-bipartite. This matrix has the following properties

* $\tilde{\Bm} \mathbf{1} = \mathbf{1}$

* It is non-symmetric, so eigenvectors  not orthogonal.

* The matrix has the spectral folding property, in anny inner produc however the all ones vector is not a generalized eigenvector (no zero DC property).

* Eigenvectors in $\Qm$ inner product, are related to eigenvectors of idea 1, similar to the relation between normalized and random walk laplacian eigenvectors.

§.§ Idea 3

Use adjacency from idea 1 and apply Sinkhorn–Knopp (SK) to make it double stochastic, then we can just do a regular normalized laplacian filterbanks. We can use the laplacian

\begin{equation}
\Id - \Cm^{-1/2}\hat{\Bm}\Cm^{-1/2}
\end{equation}

where $\Cm^{-1/2}$ is a diagonal matrix.

* Matrix is symmetrix, so orthogonal eigenvectors in standard inner product.

* Iterated filterbanks are orthogonal.

* Less interpretable, but might have decent results, given previous work that have used matrices produced by SK <cit.>.

§.§ Spectral properties of the (M,Q)-GFT

In this section we prove some spectral properties of the $(\Mm,\Qm)$-GFT that has the spectral folding property.

Eigenvalue bounds.
The generalized eigenvalues lie in $[0,2]$.
Since    $\xv^{\top}\Mm \xv \geq 0$ for all $\xv$, we also have that $\xv^{\top}\Jm \Mm \Jm \xv \geq 0$, and

\begin{equation}\label{eq_bound_xMx}
0 \leq \xv^{\top}\Mm \xv \leq \xv^{\top}\Mm \xv + \xv^{\top}\Jm \Mm \Jm \xv = 2 \xv^{\top}\Qm \xv.
\end{equation}

The last equality uses the identity $\Mm + \Jm\Mm\Jm = 2\Qm$.
Now let  $\uv$ be  a generalized  eigenvector with unit $\Qm$ norm, and eigenvalue $\lambda$, thus $\uv^{\top} \Mm \uv = \lambda \uv^{\top} \Qm \uv = \lambda$. Then using  (<ref>), we have that

\begin{equation}
0 \leq   \frac{\uv^{\top}\Mm \uv}{\uv^{\top}\Qm \uv} = \lambda  \leq 2.
\end{equation}

This implies that SGFs based on the $(\Mm,\Qm)$-GFT can be designed as continuous functions on  $[0,2]$, independently of the graph and variation operator, as it is traditionally done for the $(\Lcb, \Id)$-GFT.

Multiplicity of $\lambda=1$.
We have the  lower bound,

\begin{equation}\label{eq_lambda1_dimension}
\vert \{ i: \lambda_i=1 \} \vert \geq  \vert \vert \Ac \vert - \vert \Bc \vert   \vert.
\end{equation}

The proof of (<ref>) follows from a chain of algebraic arguments. We first use the fact that  $\Id - \Zm$ has a bipartite sparsity pattern, thus $(\Id - \Zm)^2$ is block diagonal, then

\begin{align}
\vert \{ i: \lambda_i=1 \} \vert &= \dim\ker(\Id - \Zm ) \\
&= n - \rank(\Zm_{\Ac\Bc}) - \rank(\Zm_{\Bc\Ac})\\
&= n-2\rank(\Mm_{\Ac\Bc})\\
&\geq n - 2 \min(\vert \Ac \vert, \vert \Bc \vert) \\
%&= \max(\vert \Ac \vert, \vert \Bc \vert) - \min(\vert \Ac \vert, \vert \Bc \vert) \\
&= \vert \vert \Ac \vert - \vert \Bc \vert \vert.
\end{align}

The rank of a  bipartite matrix is the sum of the ranks  of the blocks (this can be proven using the fact that the square of bipartite is block diagonal).  Then the rank of each individual term is equal to the rank of $\Mm_{\Ac\Bc}$, because $\Qm$ is non-singular. The rank of a matrix is upper bounded by the minimum between the number of rows and columns. Finally we use the fact that $n = \vert \Ac \vert + \vert \Bc \vert$.

Multiplicity of $\lambda_1$. When $\Mm \succeq 0$ is a generalized Laplacian, that is, $\Mm_{ij} \leq 0$ for all $i \neq j$, the multiplicity of $\lambda_1$ is equal to the number of connected components of the graph.
To see this note that the operation $\hat{\Mm} = \Mm - \lambda_1\Qm$ preserves eigenvectors, and shifts all eigenvalues by $-\lambda_1$. We can write $\Mm = \Qm - \Wblock$, where $\Wblock$ is a non negative matrix consisting of the block diagonal part of $\Mm$, thus $\hat{\Mm} = (1-\lambda_1)\Qm - \Wblock$.
Since $\lambda_1 \leq 1$, then $\hat\Mm$ is   a positive semidefinite generalized Laplacian matrix [71], thus $\dim\ker(\hat\Mm)$, the number of connected components, and the   multiplicity of $\lambda_1$ are equal. This property agrees with known results about the smallest eigenvalue of the combinatorial and normalized Laplacians.

§.§ Proof of <ref>

Recall that $\Qm$ from (<ref>) obeys

\begin{equation}
\Qm = \Dbip + \Lblock = \Dm - \Wblock.
\end{equation}

Using this notation, the fundamental matrix can be written

\begin{equation}
\Zm = \Id - (\Dm - \Wblock)^{-1}\Wbip.
\end{equation}

Because $\Qm^{-1}$   is non negative (see (<ref>)), and
$\Zm \mathbf{1=0}$,   the matrix $(\Dm - \Wblock)^{-1}\Wbip$ is  right stochastic.
In fact, we can write it as follows

\begin{equation}\label{eq_normalized_WQ}
(\Dm - \Wblock)^{-1}\Wbip =  \Qm^{-1}\Wbip =   \Pblock \Pbip,
\end{equation}

where $\Pbip = (\Dbip)^{-1}\Wbip$, and $\Pblock = ((\Dbip)^{-1}(\Dbip + \Lblock))^{-1}$, are right stochastic matrices. We can verify this by computing

\begin{equation}
(\Dbip)^{-1}\Wbip \mathbf{1} = (\Dbip)^{-1}\Dbip \mathbf{1} = \mathbf{1}.
\end{equation}

Note that $(\Pblock)^{-1}$ is right stochastic since

\begin{align*}
(\Dbip)^{-1}(\Dbip + \Lblock) \mathbf{1} &=    (\Dbip)^{-1}(\Dbip\mathbf{1} + \Lblock\mathbf{1})\\
&=(\Dbip)^{-1}(\Dbip\mathbf{1}) = \mathbf{1}.
\end{align*}

We have proven that $(\Pblock)^{-1} \mathbf{1}=\mathbf{1}$, and conclude by right multiplying by $\Pblock$.

§.§ Proof of Theorem <ref>

The projector  onto the eigenspace with eigenvalue $\lambda$ is

\begin{equation}
\Pm_{\lambda} = \sum_{i: \lambda_i = \lambda} \uv_i \uv_i^{\top} \Qm,
\end{equation}

It is easy to verity that $\Pm_{\lambda}^2 = \Pm_{\lambda}$, and $\Pm_{\lambda} \Pm_{\gamma} = \mathbf{0}$ when $\lambda \neq \gamma$.
Spectral graph filters can be written using the projectors,

\begin{equation}
\Hm_i = \sum_{\lambda \in \sigma(\Mm,\Qm)} h_i(\lambda) \Pm_{\lambda}, \quad \Gm_i = \sum_{\lambda \in \sigma(\Mm,\Qm)} g_i(\lambda) \Pm_{\lambda}.
\end{equation}

Then we have that

\begin{align}
&\mathbf{G}_0 \mathbf{H}_0 + \mathbf{G}_1 \mathbf{H}_1 = \sum_{\mathclap{\lambda \in \sigma(\Mm,\Qm)}} (g_0(\lambda)h_0(\lambda) + g_1(\lambda)h_1(\lambda))\Pm_{\lambda}, \\
&\mathbf{G}_0 \Jm \mathbf{H}_0 - \mathbf{G}_1\Jm \mathbf{H}_1=
\sum_{\mathclap{\lambda,\gamma \in \sigma(\Mm,\Qm)}}(g_0(\lambda) h_0(\gamma) -g_1(\lambda) h_1(\gamma) )\Pm_{\lambda}\Jm \Pm_{\gamma}.
\end{align}

<ref> implies $\Jm \Pm_{\gamma} = \Pm_{2-\gamma}\Jm$. Using  orthogonality,  the only terms that survive are those for which $\lambda = 2 -\gamma$, hence

\begin{align}
&\mathbf{G}_0 \Jm \mathbf{H}_0 - \mathbf{G}_1\Jm \mathbf{H}_1=\nonumber \\
&\sum_{\lambda \in \sigma(\Mm,\Qm)}(g_0(\lambda) h_0(2-\lambda) -g_1(\lambda) h_1(2-\lambda) )\Pm_{\lambda}\Jm.
\end{align}

The proof follows directly from these identities. Note that if (<ref>) and (<ref>) are true, then $\mathbf{G}_0 \Jm \mathbf{H}_0 - \mathbf{G}_1\Jm \mathbf{H}_1 = \mathbf{0}$, and $\mathbf{G}_0 \mathbf{H}_0 + \mathbf{G}_1 \mathbf{H}_1 = 2\sum_{\lambda \in \sigma(\Mm,\Qm)}  \Pm_{\lambda} = 2\Id$. For the converse, assume that $\mathbf{G}_0 \mathbf{H}_0 + \mathbf{G}_1 \mathbf{H}_1 + \mathbf{G}_0 \Jm \mathbf{H}_0 - \mathbf{G}_1\Jm \mathbf{H}_1 = 2\Id$. Now let $\gamma \in \sigma(\Mm,\Qm)$. Right multiplication by $\Pm_{\gamma}$ implies (<ref>), while right multiplication by $\Pm_{2-\gamma}$ leads to   (<ref>).

§.§ Proof of  <ref>

We need to show that $\Tm_a^{\top} \Qm \Tm_a = \Qm$. We will use the following properties.

* Since $\Qm$ is block diagonal, we have that $\Qm\Jm = \Jm \Qm$.

* If $\Hm$ is a spectral filter, then $\Hm^{\top}\Qm= \Qm \Hm$.  Indeed $\Hm =\Um h(\Lam)\Um^{\top}\Qm $, then

\begin{equation}
\Hm^{\top}\Qm   = \Qm \Um h(\Lam)\Um^{\top}\Qm = \Qm \Hm.
\end{equation}

Since $\Id + \Jm$ and $\Id - \Jm$ have disjoint support,

\begin{align}
\Tm_a^{\top} \Qm \Tm_a &=\frac{1}{2} \Hm_0^{\top}\Qm (\Id + \Jm)\Hm_0 +  \frac{1}{2} \Hm_1^{\top} \Qm (\Id - \Jm)\Hm_1 \\
&=\frac{1}{2}\Qm  \Hm_0 (\Id + \Jm)\Hm_0 +  \frac{1}{2} \Qm \Hm_1 (\Id - \Jm)\Hm_1.
\end{align}

Using the projectors $\Pm_{\lambda}$, we  can write

\begin{multline}
2\Tm_a^{\top} \Qm \Tm_a =\sum_{\mathclap{\lambda \in \sigma(\Mm,\Qm)}} (h_0^2(\lambda) + h_1^2(\lambda))\Qm \Pm_{\lambda} + \\
\sum_{\mathclap{\lambda \in \sigma(\Mm,\Qm)}} (h_0(\lambda)h_0(2-\lambda) - h_1(\lambda)h_1(2-\lambda))\Qm \Pm_{\lambda}\Jm.
\end{multline}

If conditions (<ref>) and (<ref>) are true, then $\Tm_a^{\top} \Qm \Tm_a = \Qm$. To show the  converse, we use   the same strategy used to prove Theorem <ref>, where we right multiply by $\Pm_{\gamma}$ and $\Pm_{2-\gamma}$.

§.§ Proof of <ref>

Let   $ \cv =[\av_0^{\top}, \dv_0^{\top}, \cdots, \dv_{L-1}^{\top}]$, and $\xv = \Tcb_s \cv$.
The reason we cannot apply (<ref>) recursively, is because (<ref>) does not hold. We can instead use Lemma <ref>, and  the inequalities

\begin{align} \label{eq_bound_norms_iteratedFB}
\alpha_{\ell} \Vert \av_{\ell} \Vert_{\Qm_{\ell}}^2 \leq  \Vert \av_{\ell} \Vert_{\Qm_{0,\ell}}^2 \leq \beta_{\ell} \Vert \av_{\ell} \Vert_{\Qm_{\ell}}^2.
\end{align}

Indeed, using $\Qm_L$ orthogonality   and (<ref>) implies

\begin{align}
\Vert \xv \Vert_{\Qm_L}^2 &= \Vert \av_{L-1} \Vert_{\Qm_{0,L-1}}^2  + \Vert \dv_{L-1} \Vert_{\Qm_{0,L-1}}^2\\
&\leq \beta_{L-1}\Vert \av_{L-1} \Vert_{\Qm_{L-1}}^2 + \Vert \dv_{L-1} \Vert_{\Qm_{0,L-1}}^2.
\end{align}

We  apply  $\Qm_{\ell}$ orthogonality and (<ref>) for each $\ell \leq L-1$ and,

\begin{equation}
\Vert \xv \Vert_{\Qm_L}^2 \leq B_0 \Vert \av_0\Vert^2_{\Qm_{0,0}} + \sum_{\ell = 0}^{L-1}B_{\ell}  \Vert \dv_{\ell}\Vert^2_{\Qm_{1,\ell}} \leq D \Vert \cv \Vert^2_{\Qm}.
\end{equation}

We can use the lower bounds in (<ref>), to get the lower bound

\begin{equation}
C \Vert \cv \Vert^2_{\Qm} \leq A_0 \Vert \av_0\Vert^2_{\Qm_{0,0}} + \sum_{\ell = 0}^{L-1}A_{\ell}  \Vert \dv_{\ell}\Vert^2_{\Qm_{1,\ell}} \leq \Vert \xv \Vert_{\Qm_L}^2.
\end{equation}

§.§ Proof of  <ref>

First, observe that   (<ref>) and (<ref>) imply that for each resolution $\ell$,  $\E[\av_\ell \dv_\ell^{\top}] = \mathbf{0}$. Linearity of the expectation implies that any linear transformation of $\av_\ell$ is also uncorrelated with $\dv_\ell$, then $\E[\av_0 \dv_\ell^{\top}] = \mathbf{0}$ and $\E[\dv_\tau \dv_\ell^{\top}] = \mathbf{0}$,

for all $\ell \geq 0$, and for all $\tau \neq \ell$. Finally,  the covariance of the detail coefficients is

\begin{equation*}
\E[\dv_\ell \dv_\ell^{\top}] = \Sigmam_{\Bc_\ell \Bc_\ell} - \Sigmam_{\Bc_\ell \Ac_\ell} \Sigmam_{\Ac_\ell \Ac_\ell}^{-1}\Sigmam_{\Ac_\ell \Bc_\ell} = ((\Mm_{\ell+1})_{\Bc_\ell \Bc_\ell})^{-1}.
\end{equation*}

We conclude by using the spectral folding property and  (<ref>), leading to $(\Mm_{\ell+1})_{\Bc_\ell \Bc_\ell} = (\Qm_{\ell+1})_{\Bc_\ell \Bc_\ell} = \Qm_{1,\ell}$.

§.§ Derivation of (<ref>)

Since $\Qm = \Pm^{1/2}(\Id - \Pm^{-1/2} \Wblock \Pm^{-1/2})\Pm^{1/2}\succ 0$, we have that $\Pm^{-1/2} \Wblock \Pm^{-1/2} \prec 1$, and   $\rho(\Ac) = \Vert \Pm^{-1/2}\Wblock \Pm^{-1/2} \Vert <1$. We have the bound $\Vert \Qm \Vert \leq \Vert \Pm \Vert(1+\rho(\Ac))$. The fact that $\rho(\Ac)<1$ also implies

\begin{equation}\label{eq_inverse_Q}
\Qm^{-1} =  \Pm^{-1/2}\left( \Id -\Pm^{-1/2}\Wblock \Pm^{-1/2}   \right)^{-1} \Pm^{-1/2}.
\end{equation}

The norm of $\Qm^{-1}$ can be bounded as $\Vert \Qm^{-1} \Vert \leq \Vert \Pm^{-1} \Vert (1-\rho(\Ac))^{-1}$. Combining the upper bounds for $\Vert \Qm \Vert$ and $\Vert \Qm^{-1}\Vert$ results in the bound for the condition number

\begin{equation}
\kappa(\Qm) \leq \kappa(\Pm) ({1+\rho(\Ac)})/({1-\rho(\Ac)}).
\end{equation}