1

Nonsubsampled Graph Filter Banks and Distributed
Implementation
Junzheng Jiang, Cheng Cheng, and Qiyu Sun

arXiv:1709.04107v1 [] 13 Sep 2017

Abstract
In this paper, we consider nonsubsampled graph filter banks (NSGFBs) to process data on a graph in a
distributed manner. Given an analysis filter bank with small bandwidth, we propose algebraic and optimization
methods of constructing synthesis filter banks such that the corresponding NSGFBs provide a perfect signal
reconstruction in the noiseless setting. Moreover, we prove that the proposed NSGFBs can control the resonance
effect in the presence of bounded noise and they can limit the influence of shot noise primarily to a small
neighborhood of its location on the graph. For an NSGFB on a graph of large size, a distributed implementation
has a significant advantage, since data processing and communication demands for the agent at each vertex depend
mainly on its neighboring topology. In this paper, we propose an iterative distributed algorithm to implement the
proposed NSGFBs. Based on NSGFBs, we also develop a distributed denoising technique which is demonstrated
to have satisfactory performance on noise suppression.

Keywords: Graph signal processing, Graph filter bank, Distributed algorithm, Noise suppression,
Random geometric graph, Laplacian matrix.
I. Introduction
Spatially distributed networks (SDNs) have an agent at each location equipped with some data processing
and communication abilities, and they have been widely used in wireless sensor networks, power grids
and many real world applications ([1]–[5]). Data collected by an SDN resides naturally on vertices of
a graph. Graph signal processing provides an innovative framework to process data on graphs. Many
concepts, such as the Fourier transform, wavelet transform and filter banks, in classical signal processing,
have been extended to graph settings in recent years. However there are still lots of fundamental problems
unexplored or not completely answered ([6]–[11]).
The wavelet transform is one of the most prominent techniques to process signals in regular domains
([12]–[14]). During the past decades, graph wavelet transforms have been introduced and some of them are
designed using the eigenvalue and eigenspace information of the graph Laplacian matrix ([15]–[19]). Graph
wavelet transform is under the same theoretical structure with graph filter banks and the corresponding
wavelet filter bank carries a down and up-sampling procedure ([7], [8], [20]–[27]). Several forms of
the down and up-sampling procedure have been defined by the partitioned graph coloring in [20], the
maximum spanning tree structure of the graph in [25], and the SVD decomposition of the graph Laplacian
matrix in [26]. A proper definition of the down and up-sampling procedure is not obvious especially when
the residing graph is of large size and complicated topological structure. This motivates us to consider a
nonsubsampled graph filter bank (NSGFB) that contains an analysis filter bank (H0 , H1 ) and a synthesis
filter bank (G0 , G1 ), see Figure 1 for its block diagram. The analysis filter bank decomposes a graph
signal into two components carrying different frequency information. The nonsubsampled structure in an
NSGFB greatly simplifies the design of analysis filter banks for spectral decomposition and synthesis
filter banks for signal reconstruction.
Jiang is with the School of Information and Communication, Guilin University of Electronic Technology, Guilin 541004, China; Cheng is
with the Department of Mathematics, Duke University, and Statistical and Applied Mathematical Sciences Institute (SAMSI), Durham, NC
27708; and Sun is with the Department of Mathematics, University of Central Florida, Orlando, Florida 32816. Emails: jzjiang@guet.edu.cn;
cheng87@math.duke.edu; qiyu.sun@ucf.edu. This work is partially supported by the National Natural Science Foundation of China (Grant
No. 61761011), SAMSI under the National Science Foundation (DMS-1638521), and the National Science Foundation (DMS-1412413).

2

Fig. 1. Block diagram of an NSGFB with analysis filter bank (H0 , H1 ) and synthesis filter bank (G0 , G1 ), where x is the input of the
NSGFB and x̃ is its output.

Filter banks can be implemented either in a centralized system or a cooperative decentralized (distributed) system. In a centralized system, a central facility receives data from agents at vertices, performs
designed data processing and sends the processed data back to agents at vertices. In a decentralized system,
the agent at each vertex has certain data processing ability to perform designed data processing, and data
collected from an agent at each vertex is shared only with neighboring vertices. Most filter banks on graphs
are designed for centralized processing, however for the implementation of filter banks on a graph of large
size, a centralized system may suffer from high computational burden and call for significant efforts to
create a data exchange network. For signal processing on an SDN or a graph of large size, a distributed
implementation provides an indispensable tool. It has been used for signal sampling and reconstruction on
an SDN in [5], graph signal inpainting in [28] and economic dispatch in power networks [29]. The reader
may refer to [5], [28]–[34] and references therein on distributed implementation of signal processing on
graphs. In a distributed implementation of the analysis and synthesis procedures of an NSGFB, signal
information on each vertex is transmitted only to neighboring vertices, which dramatically reduces the
computational cost and calls for low energy consumption. In this paper, we study NSGFBs on a cooperative
decentralized system from design to distributed implementation, and then to distributed signal denoising.
A. Main contributions
An important concept for an NSGFB is the perfect reconstruction condition, i.e., the output x̃ in Figure
1 is always the same as the input x, which can be characterized by the following matrix equation,
G0 H0 + G1 H1 = I,

(I.1)

where (H0 , H1 ) and (G0 , G1 ) are its analysis and synthesis filter bank respectively. Given an analysis
filter bank, the existence of synthesis filter banks is theoretically guaranteed so that the corresponding
NSGFB satisfies the perfect reconstruction condition (I.1) ([13], [14]). The first contribution of this paper
is that we introduce two methods to construct localized synthesis filter banks. In the first approach, the
synthesis filter bank is obtained by solving a Bezout identity for polynomials. Its bandwidth could be no
larger than the bandwidth of the analysis filter bank. In the second approach, the synthesis filter bank is
the solution of a constrained optimization problem. It does not necessarily have small bandwidth, however
it has an exponential off-diagonal decay. Consequently, the output of the corresponding localized NSGFB
suffers primarily in a small neighborhood of vertices where agents lose data processing ability and/or
communication capability.
In some real world applications of an NSGFB, the input x is the original signal xo corrupted by an
additive noise  . In addition, the subband signals z0 = H0 x and z1 = H1 x are usually processed via some
(non)linear procedure, such as hard/soft thresholding and quantization. Then the output
x̃ = G0 Ψ0 (H0 (xo +  )) + G1 Ψ1 (H1 (xo +  ))

(I.2)

of the NSGFB is no longer the original signal xo , where  is the input noise, and Ψ0 , Ψ1 are subband
processing operators. The robustness of an NSGFB is of paramount importance. For an SDN, an agent at
each vertex operates almost independently and the noise that arises at each vertex of the graph is usually
contained in some range [5]. So we may use a bounded determinstic/random noise model for NSGFBs on
a distributed system. A reasonable fidelity measure to assess the robustness of an NSGFB is the bounded

3

difference kx̃ − xo k∞ , instead of the conventional least squares error kx̃ − xo k2 , between the original
signal xo and the output signal x̃ ([5], [35], [36]). Here for 1 ≤ p ≤ ∞, `p is the space of all p-summable
sequences with norm k · kp . The second contribution of this paper is that for the NSGFB with analysis
filter bank having small bandwidth and synthesis filter banks obtained from our approaches, we establish
a quantitative estimate on the bounded difference kx̃ − xo k∞ , which is independent on the size of the
graph G. This indicates that the proposed NSGFB can control the resonance effect in the presence of
bounded additive noises.
For an NSGFB on a graph of large size, a distributed implementation may provide an indispensable
tool. The third contribution of this paper is that we propose an iterative distributed algorithm to implement
the synthesis procedure of an NSGFB rather than finding the synthesis filter bank explicitly. The keys
behind the algorithm are the decomposition (II.1) that splits the whole residing graph into a family of
overlapping subgraphs of appropriate size, and an observation that solutions of some global optimization
problem can be locally approximated by solutions of local optimization problems, when the objective
function and constraints are well localized [5]. As an application of NSGFBs, we develop a distributed
denoising technique that has satisfactory performance on noise suppression.
B. Organization
In Section II, we briefly review some fundamental concepts of graphs and introduce an overlapping graph
decomposition (II.1). In Section III, we introduce the concept of graph filters on `p , 1 ≤ p ≤ ∞, and show
that bounded filters with finite bandwidth are graph filters on `p , see Definition III.1 and Proposition III.3.
In Section IV, we discuss the analysis filter bank (H0 , H1 ) of an NSGFB which are required to have small
bandwidth, to pass/block the normalized constant signal, and to have stability on `2 , see Assumptions IV.1,
IV.2 and IV.4. We show that analysis filter banks have stability on `p for all 1 ≤ p ≤ ∞, with an estimate
on their lower and upper `p -stability bounds independent on the size of the graph, see Theorem IV.6. In
Section V, we propose an algebraic design of synthesis filter banks (G0 , G1 ) when analysis filters H0 and
H1 are polynomials of the symmetric normalized Laplacian on the graph, see Theorem V.1. In Section VI,
we consider the construction of synthesis filter banks (G0 , G1 ) by solving the constrained optimization
problem (VI.1) and (VI.2) with the objective function consisting of Frobenius norms of G0 and G1 , see
(VI.4) and Theorem VI.1. In Section VII, we propose an exponentially convergent iterative algorithm
(VII.9) and (VII.10) to implement the synthesis procedure, where each iteration can be implemented in
a distributed manner, see Theorem VII.2 and Algorithm VII.1. In Section VIII, we create a distributed
denoising technique associated with spline NSGFBs, see Figure 5, and demonstrate its performance for
signal denoising on graphs of large size. All proofs are collected in the appendices.
C. Notation
We use the common convention of representing matrices and vectors with bold letters and scalars
with normal letters. For a matrix A, denote its transpose, trace, Frobenius norm and operator norm on
`p , 1 ≤ p ≤ ∞, by AT , tr(A), kAkF and kAkBp respectively. For a graph G, denote its adjacency matrix
and degree matrix by AG and DG respectively, and define its Laplacian matrix by LG := DG − AG and
−1/2
−1/2
its symmetric normalized Laplacian matrix by Lsym
:= DG LG DG . For a scalar t, let sgn(t), btc
G
and t+ be its sign, integral part and positive part respectively, and t be the vector of appropriate size
with all entries taking value t. For a set F , denote its cardinality and indicator function by #F and χF
respectively.
II. P RELIMINARIES ON GRAPHS
Let G := (V, E) be a graph, where V = {1, 2, · · · , N } is the set of vertices and E is the set of edges
([5], [6]). For the distributed implementation of an NSGFB, we require that the residing graph G has
certain global features.

4

Assumption II.1. Throughout the paper, we consider simple graphs G, i.e., they are undirected and
unweighted, and they do not contain self-loops and multiple edges.
Take r ≥ 0. For a graph G = (V, E) satisfying Assumption
II.1, we define the r-neighborhood and

the r-neighboring subgraph of i ∈ V by B(i, r) := j ∈ V : ρ(i, j) ≤ r and Gi,r := B(i, r), E(i, r)
respectively, where E(i, r) contains all edges of the graph G with endpoints in B(i, r), and ρ(i, j) is the
geodesic distance between vertices i and j in V . Then for r ≥ 1, we can decompose the graph G into a
family of overlapping subgraphs Gi,r , i ∈ V , of diameters at most 2r,
G = ∪i∈V Gi,r .

(II.1)

For distributed implementation for an NSGFB, we presume that numbers of vertices in the r-neighborhood
of any vertex are dominated by a polynomial about r.
Assumption II.2. Throughout the paper, we consider graphs G with the counting measure µ having
polynomial growth, i.e., there exist positive constants D1 (G) and d such that
µ(B(i, r)) ≤ D1 (G)(r + 1)d

(II.2)

for all i ∈ V and r ≥ 0, where µ(F ) := #F for all F ⊂ V .
The minimal constants d and D1 (G) in (II.2) are called as Beurling dimension and density of the graph
G respectively [5].
The decomposition (II.1) plays a crucial role in the proposed distributed implementation for an NSGFB,
and the selection of the radius parameter r in (II.1) depends on Beurling dimension d and density D1 (G)
of the graph G, see Theorem VII.2. Accordingly, we expect that the Beurling dimension d and density
D1 (G) of the graph G are much smaller than (or even independent on) the size of the graph, which
implies that the graph G should be sparse. Shown in Figure 2 are two representative graphs that satisfy
Assumptions II.1 and II.2:
• The Minnesota traffic graph with 2642 vertices, where each vertex represents a spatial location in
the state of Minnesota equipped with a traffic monitoring sensor and each edge denotes a direct
communication link between monitoring sensors ([20], [21]).
• The random geometric graph RGGN with N vertices randomly deployed in the√
region [0, 1]2 and an
edge existing between two vertices if their physical distance is not larger than 2N −1/2 ([26], [37],
[38]).

Fig. 2. Plotted on the left is the Minnesota traffic graph that has Beurling dimension 2 and density 2.1378. On the right is a random
geometric graph with N = 4096, which has Beurling dimension 2 and density 3.0775.

5

III. G RAPH SIGNAL AND FILTERING
Let G = (V, E) satisfy Assumptions II.1 and II.2. A signal x residing on the graph G is a vector
(xi )i∈V , where xi refers to the signal value at vertex i ∈ V . In SDNs and many real world applications,
data collected belongs to some sequence space `p , 1 ≤ p ≤ ∞ ([5], [35], [36]).
A filter A on the graph G is a linear transformation from one signal x on G to another signal y = Ax
on G, which is usually represented by a matrix A = (a(i, j))i,j∈V . A filter A is expected to map a signal
with finite energy to another signal with finite energy and a bounded signal to another bounded signal. A
quantitative description of the above filtering procedure is
kAxkp ≤ Ckxkp for all x ∈ `p ,

(III.1)

where 1 ≤ p ≤ ∞ and C is a positive constant.
Definition III.1. Let 1 ≤ p ≤ ∞. We say that A is a graph filter on `p if (III.1) is satisfied, and we call
the minimal constant C for (III.1) to hold, denoted by kAkBp , the filter bound on `p .
In some
P practical applications ([8], [20], [21], [22], [26], [39]), a graph filter A is a polynomial
P (t) := Ll=0 pl tl of the symmetric normalized Laplacian Lsym
on G, i.e.,
G
A=

P (Lsym
G )

= p0 I +

L
X

l
pl (Lsym
G ).

(III.2)

l=1

Let 0 ≤ λ1 ≤ λ2 ≤ · · · ≤ λN ≤ 2 be eigenvalues of the symmetric normalized Laplacian Lsym
and write
G
Lsym
= UT Λ U,
G

(III.3)

where UT = [u1 , . . . , uN ] is an orthogonal matrix and Λ = diag(λ1 , . . . , λN ) is a diagonal matrix. Then
Λ)U,
A = UT P (Λ

(III.4)

and the filter bound kAkB2 can be evaluated explicitly,
kAkB2 = sup |P (λn )| ≤ sup |P (t)|.

(III.5)

0≤t≤2

1≤n≤N

To estimate kAkBp , p 6= 2, of a graph filter A = (a(i, j))i,j∈V , we define the bound of A by
kAk∞ = sup |a(i, j)|.

(III.6)

i,j∈V

A graph filter A on `p , 1 ≤ p ≤ ∞, has bounded entries and
kAk∞ ≤ sup kAej kp ≤ kAkBp , 1 ≤ p ≤ ∞,

(III.7)

j∈V

where the last inequality is obtained from (III.1) by replacing x by the standard unit vector ej with j-th
component taking value one while all others are zero.
For the distributed implementation for an NSGFB, bounded filters with finite bandwidth will be used
as analysis filters, see Assumption IV.1.
Definition III.2. The bandwidth σ := σ(A) of a graph filter A = (a(i, j))i,j∈V is the minimal nonnegative
integer such that a(i, j) = 0 for all i, j ∈ V with ρ(i, j) > σ. For a filter pair (A, B), we define its
bandwidth σ := σ(A, B) by max(σ(A), σ(B)).
In the following proposition, we show that a bounded filter with finite bandwidth is a graph filter on
`p , 1 ≤ p ≤ ∞, with filter bound dominated by some constant, independent of the size of the graph G.

6

Proposition III.3. Let 1 ≤ p ≤ ∞, and G be a graph satisfying Assumptions II.1 and II.2. Then for any
bounded graph filter A with bandwidth σ, we have
kAk∞ ≤ kAkBp ≤ D1 (G)(σ + 1)d kAk∞ ,

(III.8)

where d and D1 (G) are the Beurling dimension and density of the graph G respectively.
spln
For n ≥ 1, we define spline filters Hspln
0,n and H1,n of order n by
n
1

1 sym n
spln
sym
spln
and H1,n =
L
,
H0,n = I − LG
2
2 G

(III.9)

spln
see [39] in circulant graph setting. Spline filters Hspln
0,n and H1,n , n ≥ 1, have bandwidth n and their filter
bounds on `2 dominated by one, i.e.,

kHspln
l,n kB2 ≤ 1, l = 0, 1.

(III.10)

For p 6= 2, we obtain from Proposition III.3 that
spln n
spln
d
n
kHspln
l,n kBp ≤ kHl,1 kBp ≤ 2 D1 (G)kHl,1 k∞ )

= 2d−1 D1 (G))n , l = 0, 1,

(III.11)

where d and D1 (G) are the Beurling dimension and density of the graph G respectively. Therefore our
spln
p
representative spline filters Hspln
0,n and H1,n , n ≥ 1, are graph filters on ` , 1 ≤ p ≤ ∞, with filter bounds
dominated by some constants independent of the size of the graph G.
IV. A NALYSIS FILTER BANKS
The analysis filter bank decomposes the input signal on a graph into two components carrying frequency
information. In this section, we design the analysis filter bank (H0 , H1 ) of an NSGFB to have small
bandwidth, to pass/block the normalized constant signal, and to have stability on `2 , see Assumptions
IV.1, IV.2 and IV.4. In this section, we also show that analysis filter banks have stability on `p for all
1 ≤ p ≤ ∞, with an estimate on their lower and upper `p -stability bounds independent of the size of the
graph, see Theorem IV.6.
Let G = (V, E) be a graph satisfying Assumptions II.1 and II.2, and (H0 , H1 ) be the analysis filter
bank of an NSGFB. For the distributed implementation of an NSGFB, we make the following assumption
for its analysis filter bank (H0 , H1 ).
Assumption IV.1. The analysis filter bank (H0 , H1 ) has bandwidth σ ≥ 1.
Given an input graph signal x = (xi )i∈V , outputs of analysis procedure are
z0 = H0 x and z1 = H1 x.

(IV.1)

Write zl = (zl (i))i∈V and Hl = (hl (i, j))i,j∈V , l = 0, 1. Then it follows from (IV.1) and Assumption IV.1
that component values of the outputs z0 and z1 at each vertex k ∈ V are weighted sums of values of the
input x in a σ-neighborhood of k,
X
zl (k) =
hl (k, i)x(i), i ∈ V.
(IV.2)
ρ(i,k)≤σ

Thus the analysis procedure of an NSGFB can be implemented in a distributed manner.
To apply an NSGFB to some real world applications, such as noise suppression and abnormal phenomenon detection, its analysis filter bank should constitute certain spectral decomposition ([20], [21],
[23], [45]). Throughout the paper, we also make the following assumption for the analysis filter bank
(H0 , H1 ).

7

1/2

Assumption IV.2. The filter H0 passes the normalized constant signal DG 1, and the filter H1 blocks
1/2
the normalized constant signal DG 1, i.e.,
1/2

1/2

1/2

H0 DG 1 = DG 1 and H1 DG 1 = 0.

(IV.3)

The frequency partition of an analysis filter bank on an arbitrary graph G is not as obvious as that in
classical setting. For the case that
sym
H0 = P0 (Lsym
G ) and H1 = P1 (LG )

(IV.4)

for some polynomials P0 and P1 , one may verify that Assumption IV.2 is satisfied if and only if
P0 (0) = 1 and P1 (0) = 0.

(IV.5)

1/2

The above equivalence follows from the fact that DG 1 is an eigenvector of the symmetric normalized
Laplacian Lsym
associated with eigenvalue zero.
G
spln
n
The spline filter banks (Hspln
0,n , H1,n ), n ≥ 1, are of the form (IV.4) with P0 (t) = (1 − t/2) and
P1 (t) = (t/2)n , and they satisfy Assumption IV.2 by (IV.5), i.e.,
1/2

1/2

1/2

spln
Hspln
0,n DG 1 = DG 1 and H1,n DG 1 = 0.

(IV.6)

Spline filter banks in the circulant graph setting are known in [39] as graph-spline wavelet transform.
spln
Shown in Figure 3 is local smoothing/blocking phenomenon of the spline filter bank (Hspln
0,2 , H1,2 ) to a
blockwise constant signal on the Minnesota traffic graph and a blockwise smooth signal on the random
geometric graph RGG4096 in Figure 2. It is observed that the lowpass filtered signal is very close to the
original signal except near the boundary between different blocks, and that the highpass filtered signal
essentially vanishes except around the region where the original signal exhibits sharp local variation.

Fig. 3. Plotted on the top (resp. at the bottom), from left to right, are the original signal x on the Minnesota traffic graph (resp. on the
spln
random geometric graph RGG4096 in Figure 2), the lowpass filtered signal Hspln
0,2 x and the highpass filtered signal H1,2 x. The signal x on
the top is a blockwise constant function that has only two values ±1 on three blocks with one block only containing a vertex ([20], [21]),
and the signal x at the bottom is a blockwise polynomial consisting of four strips and imposing the polynomial 0.5 − 2cx on the first and
third diagonal strips and 0.5 + c2x + c2y on the second and fourth strips respectively, where (cx , cy ) are the coordinates of vertices ([26]).

Robustness is a fundamental requirement in the context of filter bank to control the signal dynamic
range and to regulate the input noise. For the robustness of an NSGFB on `p , 1 ≤ p ≤ ∞, we introduce
stability of a graph filter pair on `p .
Definition IV.3. Let 1 ≤ p ≤ ∞. We say that (H0 , H1 ) is stable on `p if there are two positive constants
Cp and Dp such that
1/p
Cp kxkp ≤ kH0 xkpp + kH1 xkpp
≤ Dp kxkp
(IV.7)

8

hold for all x ∈ `p if 1 ≤ p < ∞, and

C∞ kxk∞ ≤ max(kH0 xk∞ , kH1 xk∞ ≤ D∞ kxk∞

(IV.8)

hold for all x ∈ `∞ if p = ∞. The optimal constants Cp and Dp for the inequalities in (IV.7) and (IV.8) to
hold are called as lower and upper stability bounds of the graph filter bank (H0 , H1 ) on `p respectively.
Given an NSGFB with the analysis filter bank (H0 , H1 ) and synthesis filter bank (G0 , G1 ) such that
the perfect reconstruction condition (I.1) holds, we have that
(kG0 k2B2 + kG1 k2B2 )−1 kxk22 ≤ kH0 xk22 + kH1 xk22
≤ (kH0 k2B2 + kH1 k2B2 )kxk22
hold for all x ∈ `2 . So throughout the paper, we assume that the analysis procedure is stable on `2 .
Assumption IV.4. The analysis filter bank (H0 , H1 ) is stable on `2 .
For any x ∈ `2 , direct calculation leads to

kH0 xk22 + kH1 xk22 = xT (HT0 H0 + HT1 H1 x.

(IV.9)

Thus we have the following characterization to Assumption IV.4.
Proposition IV.5. Let G satisfy Assumptions II.1 and II.2. Then (H0 , H1 ) satisfies Assumption IV.4 if
and only if H := HT0 H0 + HT1 H1 is positive definite. Moreover, the optimal constants C2 and D2 for
(IV.7) to hold can be evaluated by
C22 = (kH−1 kB2 )−1 and D22 = kHkB2 .
For graph filters H0 and H1 of the form (IV.4), we obtain from (III.3) that

Λ))2 + (P1 (Λ
Λ))2 U.
HT0 H0 + HT1 H1 = UT (P0 (Λ

(IV.10)

(IV.11)

Hence we can evaluate the optimal constants C2 and D2 for (IV.7) to hold explicitly:
inf (P0 (λm ))2 + (P1 (λm ))2

1≤m≤N

= inf kH0 xk22 + kH1 xk22 ≤ sup kH0 xk22 + kH1 xk22
kxk2 =1

2

kxk2 =1
2

= sup (P0 (λm )) + (Q0 (λm )) .

(IV.12)

1≤m≤N

Set Rn (t) = (1 − t/2)2n + (t/2)2n , n ≥ 1. Then
inf Rn (t) = 2−2n+1 and sup Rn (t) = 1.

0≤t≤2

(IV.13)

0≤t≤2

Taking P0 (t) = (1 − t/2)n and P1 (t) = (t/2)n in (IV.12) and applying (IV.13), we get

spln
spln T
spln
T
2−2n+1 kxk22 ≤ xT (Hspln
)
H
+
(H
)
H
0,n
0,n
1,n
1,n x
spln
2
2
2
= kHspln
0,n xk2 + kH1,n xk2 ≤ kxk2

(IV.14)

spln
for all x ∈ `2 . Therefore spline filter banks (Hspln
0,n , H1,n ) of order n ≥ 1 satisfy Assumption IV.4 with
lower bound 2−n+1/2 and upper bound 1 by (IV.12) and (IV.14).
Filters in a stable filter bank on `p are graph filters on `p , 1 ≤ p ≤ ∞. In the following theorem, we
show that analysis filter banks are stable on `p , 1 ≤ p ≤ ∞, with quantitative estimates on their lower
and upper stability bounds by some constants independent of the size of the graph.

9

Theorem IV.6. Let G be a graph satisfying Assumptions II.1 and II.2, H0 and H1 have bandwidth σ ≥ 1,
and set H := HT0 H0 + HT1 H1 . If (H0 , H1 ) is stable on `2 , then it is stable on `p for all 1 ≤ p ≤ ∞.
Moreover, we have the following estimates for its lower and upper stability bounds Cp and Dp :
1/2

kHkB2
Cp ≥
d+1
d!2 (D1 (G))2 (σ + 1)2d κd+2
and

1/2

Dp ≤ 2D1 (G)(σ + 1)d kHkB2 ,

(IV.15)

(IV.16)

where d and D1 (G) are the Beurling dimension and density of the graph G respectively, and
κ = kH−1 kB2 kHkB2

(IV.17)

is the condition number of the matrix H.
spln
p
Combining (IV.14) and Theorem IV.6, the spline filter banks (Hspln
0,n , H1,n ), n ≥ 1, are stable on ` , 1 ≤
p ≤ ∞, and their lower and upper stability bounds Cp and Dp satisfy

22(d+2)n−1 d!(n

1
≤ Cp ≤ Dp
+ 1)2d (D1 (G))2
≤ 2D1 (G)(n + 1)d+1 .

We finish this section with a remark on stability bounds of a graph filter bank on the space `2 and on
the spaces `p , p 6= 2.
Remark IV.7. For a finite graph G = (V, E), a stable filter bank (H0 , H1 ) on `2 is also stable on
`p , 1 ≤ p ≤ ∞, and the lower stability bounds C2 and Cp satisfy
N −|1/p−1/2| ≤

C2
≤ N |1/p−1/2| ,
Cp

(IV.18)

where N = #V is the size of the graph G. The above estimation is unfavorable when the graph G has
large size, however it cannot be improved if there is no restriction on the filter bank (H0 , H1 ). As our
analysis filter bank (H0 , H1 ) has small bandwidth σ, we obtain the following estimate independent of the
size N of the graph G from Proposition IV.5 and Theorem IV.6,
1
C2
≤
≤ d!2d+1 (D1 (G))2 (σ + 1)2d κd+3/2 ,
d
1/2
2D1 (G)(σ + 1) κ
Cp

(IV.19)

where κ is given in (IV.17). The reader may refer to [5] and [40]–[44] for historical remarks and various
estimates on the ratio between stability bounds on `p and `q , 1 ≤ p, q ≤ ∞, for matrices with certain
off-diagonal decay.
V. S YNTHESIS FILTER BANKS AND B EZOUT IDENTITY
Let G = (V, E) be a graph satisfying Assumptions II.1 and II.2, and (H0 , H1 ) be a graph filter bank
satisfying Assumptions IV.1, IV.2 and IV.4. In this section, we propose an algebraic method to construct
graph filters G0 and G1 so that the NSGFB with the analysis filter bank (H0 , H1 ) and synthesis filter
bank (G0 , G1 ) satisfies the perfect reconstruction condition (I.1) and the bandwidth of synthesis filter bank
(G0 , G1 ) is no larger than the bandwidth of the analysis filter bank (H0 , H1 ). The proposed approach
applies for filter banks (H0 , H1 ) being polynomials of the symmetric normalized Laplacian on the graph
G, i.e.,
sym
H0 = P0 (Lsym
(V.1)
G ) and H1 = P1 (LG )
for some polynomials P0 and P1 .

10

Theorem V.1. Let G be a graph satisfying Assumptions II.1 and II.2, (H0 , H1 ) be a graph filter bank
satisfying Assumptions IV.1, IV.2 and IV.4 and being of the form (V.1), and let 0 ≤ λ1 ≤ λ2 ≤ · · · ≤
λN ≤ 2 be eigenvalues of the symmetric normalized Laplacian Lsym
G . If polynomials Q0 and Q1 satisfy
P0 (λm )Q0 (λm ) + P1 (λm )Q1 (λm ) = 1, 1 ≤ m ≤ N,

(V.2)

then the NSGFB with the analysis filter bank (H0 , H1 ) and synthesis filter bank (G0 , G1 ) satisfies the
perfect reconstruction condition (I.1), where
sym
G0 = Q0 (Lsym
G ) and G1 = Q1 (LG ).
1/2

(V.3)
1/2

1/2

The filter G0 in (V.3) passes the normalized constant signal DG 1, since G0 DG 1 = Q0 (0)DG 1 =
1/2
DG 1, where the last equation follows from (IV.5) and (V.2). However, the filter G1 in (V.3) may not
1/2
1/2
block the normalized constant signal, as G1 DG 1 = Q1 (0)DG 1 is not necessarily a zero signal. In this
case, we can construct a new synthesis filter bank by lifting,
e 0 = G0 + Q1 (0)H1 and G
e 1 = G1 − Q1 (0)H0 ,
G

(V.4)

e 0 D1/2 1 = D1/2 1 and G
e 1 D1/2 1 = 0.
which satisfies G
G
G
G
A strong version of (V.2) is the Bezout identity
P0 (z)Q0 (z) + P1 (z)Q1 (z) = 1, z ∈ C

(V.5)

for polynomials P0 , P1 , Q0 and Q1 . In the circulant graph setting, the above approach of constructing
synthesis filter banks via solving Bezout identity (V.5) was discussed in [39]. Comparing with the Bezout
identity (V.2) on the eigenvalue set of the symmetric normalized Laplacian matrix Lsym
G , the advantage
of the approach (V.5) provides a tool to design synthesis filter banks without a priori knowledge of
global topology of the residing graph and then it simplifies the design of synthesis filter banks for signal
reconstruction. It is well known that the Bezout identity (V.5) is solvable if and only if polynomials P0
and P1 have no common root. Moreover, there is a unique solution pair (Q0,B , Q1,B ) to the Bezout identity
(V.5) such that Q0,B (0) = 1, Q1,B (0) = 0 and the degree of Q0,B (resp. Q1,B ) is no larger than the degree
of P1 (resp. P0 ). Define
sym
G0,B = Q0,B (Lsym
(V.6)
G ) and G1,B = Q1,B (LG ).
Then the bandwidth of the synthesis filter bank (G0,B , G1,B ) is no larger than bandwidth of the analysis
filter bank (H0 , H1 ). Moreover, for any synthesis filter bank (G0 , G1 ) there exists a polynomial R such
that
sym
G0 = G0,B + R(Lsym
(V.7)
G )H1 and G1 = G1,B − R(LG )H0
satisfies the perfect reconstruction condition (I.1). We remark that the above polynomials R could be
appropriately chosen for real world applications of an NSFGB.
of order n ≥ 1 by
Following (V.6), we define synthesis spline filters GB,spln
and GB,spln
0,n
1,n
B,spln
B,spln
(Lsym
= Q1,n
(Lsym
GB,spln
= QB,spln
0,n
0,n
G ) and G1,n
G ),

where
QB,spln
(t)
0,n


n−1 
X
2n − 1 
t n−1−l  t l
=
1−
2
2
l
l=0

 
2n − 1
t n
+
n−1
2

(V.8)

11

and
QB,spln
(t)
1,n

=


n−1 
X
2n − 1  t n−1−l 
l
2


2n − 1
t n
.
−
1−
2
n−1
l=0

1−

t l
2

1/2

blocks the
passes the normalized constant signal DG 1, the filter GB,spln
For n ≥ 1, the filter GB,spln
1,n
0,n
1/2
spln
normalized constant signal DG 1, and the NSGFB with the analysis spline bank (H0,n , Hspln
1,n ) and
B,spln
B,spln
synthesis filter bank (G0,n , G1,n ) satisfies the perfect reconstruction condition (I.1). The first two
B,spln
results follow from QB,spln
(0) = 1 and Q1,n
(0) = 0, while the perfect reconstruction conclusion holds
0,n
since

 t n
t n B,spln
B,spln
1−
Q0,n (t) +
Q1,n
(t)
2
2

n−1 
X
2n − 1
=
(1 − u)2n−1−l ul
l
l=0

n−1 
X
2n − 1
+
(1 − u)l u2n−1−l
l
l=0
= ((1 − u) + u)2n−1 = 1,
where u = t/2.
In real world applications of an NSGFB such as the proposed distributed denoising in Section VIII,
the subband signals z0 and z1 in (IV.1) are processed via some (non)linear procedure, such as hard/soft
thresholding and quantization. In this case, the reconstructed signal x̃ is not necessarily the same as the
original signal x. In the following theorem, we show that the difference is mainly dominated by the error
caused by the subband processing.
Proposition V.2. Let the graph G, the analysis filter bank (H0 , H1 ) and the synthesis filter bank (G0 , G1 )
be as in Theorem V.1. Assume that the error caused by the subband processing Ψl on subband signals
zl = Hl x, l = 0, 1, is dominated by  for any input signal x ∈ `p , i.e.,
kzl − Ψl (zl )kp ≤ , l = 0, 1,

(V.9)

where  ≥ 0 and 1 ≤ p ≤ ∞. For the input signal x ∈ `p , the reconstructed signal x̃ = G0 Ψ0 (z0 ) +
G1 Ψ1 (z1 ) via the corresponding NSGFB belongs to `p as well. Moreover
kx̃ − xkp ≤ D1 (G)(σ̃ + 1)d (kG0 k∞ + |G1 k∞ ),

(V.10)

where d and D1 (G) are the Beurling dimension and density of the graph G respectively, and σ̃ is the
bandwidth of the synthesis filter bank (G0 , G1 ).
We finish this section with a distributed implementation of the NSGFB with analysis/synthesis filter
banks selected in Theorem V.1. Write Gl = (gl (i, j))i,j∈V , l = 0, 1. As the synthesis filters G0 and G1
have finite bandwidth σ̃, the synthesis procedure can be implemented in a distributed manner,
X
x̃k =
(g0 (k, j)z̃0 (j) + g1 (k, j)z̃1 (j)), k ∈ V,
(V.11)
ρ(j,k)≤σ̃

where x̃ = (x̃i )i∈V is the reconstructed signal and Ψl (zl ) = (z̃l (i))i∈V , l = 0, 1, are outputs of subband
processing. Hence values of the reconstructed signals x̃ at each vertex k ∈ V are weighted sums of
values of the subband processed outputs Ψ0 (z0 ) and Ψ1 (z1 ) in a σ̃-neighborhood of k ∈ V , cf. (IV.2) for
distributed implementation of the analysis procedure.

12

Our representative subband processing procedures Ψ are hard(soft) thresholding and uniform quantization. For those cases, the subband processing Ψ is of the form Ψ(z) = (ψ(zi ))i∈V for z = (zi )i∈V , where
ψ is the hard(soft) thresholding and uniform quantization function. Thus the subband processing can be
implemented in a distributed manner and the error resulted are bounded (i.e., (V.9) holds for p = ∞) by
the hard(soft) thresholding and quantization level. This together with (IV.2) and (V.11) implies that the
NSGFB with analysis/synthesis filter banks in Theorem V.1 can be implemented in a distributed manner
too, provided that the subband processing can be.
VI. S YNTHESIS FILTER BANK AND OPTIMIZATION
Let G = (V, E) be a graph satisfying Assumptions II.1 and II.2, and (H0 , H1 ) be a graph filter bank
satisfying Assumptions IV.1, IV.2 and IV.4. In this section, we consider the construction of synthesis filter
banks (G0 , G1 ) of an NSGFB by solving the minimization problem:
minimize kG0 k2F + kG1 k2F

(VI.1)

G0 ,G1

subject to the perfect reconstruction condition
G0 H0 + G1 H1 = I.

(VI.2)

Define the Lagrange function L of the constrained optimization problem (VI.1) and (VI.2) by
L(G0 , G1 , Θ ) = kG0 k2F + kG1 k2F

ΘT .
− tr (G0 H0 + G1 H1 − I)Θ
By direct calculation, we have




∂L
= 2G0 − Θ HT0
∂G0
∂L
= 2G1 − Θ HT1
∂G1
∂L
= G0 H0 + G1 H1
∂Θ

(VI.3)
− I.

Set H = HT0 H0 + HT1 H1 . Solving

∂L
∂L
∂L
=
=
=0
Θ
∂G0
∂G1
∂Θ
leads to the unique solution of the constrained optimization problem (VI.1) and (VI.2),
G0,L = H−1 HT0 and G1,L = H−1 HT1 .

(VI.4)

The synthesis filter bank (G0,L , G1,L ) in (VI.4) satisfies
G0,L H0 + G1,L H1 = I,
1/2

and the filter G0,L passes the normalized constant signal DG 1, since
1/2

1/2

1/2

G0,L DG 1 = H−1 (HT0 H0 + HT1 H1 )DG 1 = DG 1.
1/2

We remark that G1,L may not block the normalized constant signal DG 1.
For the case that H is a diagonal matrix, the synthesis filter bank (G0,L , G1,L ) in (VI.4) has the same
bandwidth as the analysis filter bank (H0 , H1 ), and

kH−1 kB2 kHl k∞ if ρ(i, j) ≤ σ
|gl,L (i, j)| ≤
(VI.5)
0
otherwise,
where Gl,L := gl,L (i, j))i,j∈V , l = 0, 1.
Let κ be the condition number of the matrix H in (IV.17). It is well known that κ > 1 when H is not
a diagonal matrix. For κ > 1, the synthesis filter bank (G0,L , G1,L ) in (VI.4) does not necessarily have a
small bandwidth, however it always has an exponential off-diagonal decay.

13

Theorem VI.1. Let G = (V, E) be a graph satisfying Assumptions II.1 and II.2, (H0 , H1 ) be a graph
filter bank satisfying Assumptions IV.1, IV.2 and IV.4, κ be the condition number of the matrix H :=
HT0 H0 + HT1 H1 , and let Gl,L := gl,L (i, j))i,j∈V , l = 0, 1, be as in (VI.4). Assume that κ > 1, then
|gl,L (i, j)| ≤ D1 (G)(σ + 1)d (1 − 1/κ)−1/2


θ
ρ(i, j)
(VI.6)
×kH−1 kB2 kHl k∞ exp −
2σ
hold for all i, j ∈ V and l = 0, 1, where θ = ln(κ/(κ − 1)), σ ≥ 1 is the bandwidth of the analysis filter
bank (H0 , H1 ), and d and D1 (G) are the Beurling dimension and density of the graph G respectively.
Remark VI.2. Agents located at some vertices may lose data processing ability and/or communication
capability. In that case, outputs of the analysis procedure of an NSGFB can be considered as being
corrupted by shot noise. The exponential off-diagonal decay property in Theorem VI.1 implies that the
reconstructed signal suffers mainly in their neighborhood of limited size. This means that the proposed
NSGFB can limit the influence of shot noise essentially to their small neighborhoods on the graph.
Remark VI.3. By the exponential off-diagonal decay property in Theorem VI.1, the synthesis filters
(G0,L , G1,L ) are filters on `p , 1 ≤ p ≤ ∞,
kGl,L kBp ≤ d!2d (D1 (G))2 (σ + 1)2d κd+1 (1 − 1/κ)−1/2
×kH−1 kB2 kHl k∞ , l = 0, 1.

(VI.7)

The above conclusion with p = ∞ indicates that the NSGFB does not have a resonance effect.
Applying similar argument used in the proof of Proposition V.2, we have
Corollary VI.4. Let G, (H0 , H1 ), (G0,L , G1,L ) be as in Theorem VI.1, and let p, Ψ0 , Ψ1 ,  be as in
Proposition V.2. Assume that the input signal x of the corresponding NSGFB belongs to `p , then the
reconstructed signal x̃ = G0,L Ψ0 (H0 x) + G1,L Ψ1 (H1 x) via the NSGFB belongs to `p and
kx̃ − xkp ≤ d!2d (D1 (G))2 (σ + 1)2d κd+1 kH−1 kB2
×(1 − 1/κ)−1/2 kH0 k∞ + kH1 k∞ ).
(VI.8)
Solving the constrainted optimization program (VI.1) and (VI.2) associated with the analysis spline
spln
L,spln
L,spln
filter banks (Hspln
0,n , H1,n ), we obtain the synthesis spline filter bank (G0,n , G1,n ), n ≥ 1, where


 −1 spln
L,spln
spln 2
spln 2
Gl,n = H0,n + H1,n
Hl,n , l = 0, 1.
(VI.9)
L,spln
The synthesis spline filters GL,spln
and G1,n
, n ≥ 1, have full bandwidth, however they have exponential
0,n
L,spln
(i,
j))i,j∈V , l = 0, 1. By (VI.6) and Theorem VI.1, we obtain
off-diagonal decay. Write GL,spln
=
(g
l,n
l,n
that
L,spln
|gl,n
(i, j)| ≤ 23n−3/2 (22n−1 − 1)−1/2 (n + 1)d D1 (G)

 ln(22n−1 /(22n−1 − 1))
ρ(i, j)
× exp −
2n
hold for all i, j ∈ V and l = 0, 1.
λ) to describe frequency response of a filter A = P (Lsym
By (III.4), we may use P (λ
G ) of the form
(III.2), where the vector λ = (λ1 , . . . , λN ) is composed of eigenvalues 0 ≤ λ1 ≤ λ2 ≤ · · · ≤ λN ≤ 2
of the symmetric normalized Laplacian Lsym
G . Shown in Figure 4 are frequency responses of the analysis
spln
spln
spline filter banks (H0,n , H1,n ) of order n, the synthesis spline filter banks (GB,spln
, GB,spln
) in (V.8),
0,n
1,n
L,spln
L,spln
and the synthesis spline filter banks (G0,n , G1,n ) just constructed, where n = 1, 2. It is observed
spln
that the frequency responses of analysis spline filter banks (Hspln
0,n , H1,n ) and synthesis spline filter
L,spln
L,spln
banks (G0,n , G1,n ) have certain complementary property, while the synthesis spline filter banks
(GB,spln
, GB,spln
) constructed via solving a Bezout identity do not.
0,n
1,n

14

Fig. 4. Plotted on the top (resp. at the bottom) are the frequency responses of analysis/synthesis spline filters of order n on the Minnesota
traffic graph (resp. on the random geoemetric graph RGG4096 in Figure 2), where n = 1 for the left figure and n = 2 for the right figure.

VII. I TERATIVE DISTRIBUTED ALGORITHM FOR SYNTHESIS PROCEDURE
For the NSGFB with synthesis filter banks in Theorem V.1, the distributed implementation of the
corresponding synthesis procedure has been discussed in (V.11).
For the NSGFB with the analysis filter bank (H0 , H1 ) and synthesis filter bank (G0,L , G1,L ) obtained
from solving the constrained optimization problem (VI.1) and (VI.2), output x̃ of the synthesis procedure
is
x̃ = G0,L z̃0 + G1,L z̃1 ,
(VII.1)
where z̃0 and z̃1 be outputs of subband processing. As filters G0,L and G1,L may have full bandwidth, it is
infeasible to evaluate G0,L z̃0 and G1,L z̃1 directly in a distributed manner. In this paper, we do not intend
to find synthesis filters G0,L and G1,L explicitly, instead we propose an iterative distributed algorithm to
implement the synthesis procedure (VII.1).
The proposed iterative distributed algorithm is based on two pivoting observations. The first observation
is that the output signal x̃ in (VII.1) is the unique solution of the following global least squares problem:
min kH0 x − z̃0 k22 + kH1 x − z̃1 k22 ,
x

(VII.2)

which follows from (VI.4). To solve the global optimization problem (VII.2) in a distributed manner, we
introduce a family of local least squares problems,
2
2r
2
min kH0 χ2r
k x − z̃0 k2 + kH1 χk x − z̃1 k2 , k ∈ V,
x

(VII.3)

where χrk , k ∈ V , are truncation operators defined by

χrk : (x(i))i∈V 7−→ x(i)χB(k,r) (i) i∈V ,

(VII.4)

and r ≥ 1 is a radius parameter to be determined later [5]. One may verify that given any k ∈ V , the
unique solution of the local optimization problem (VII.3) is given by

2r
2r −1 2r
vk,r = χ2r
χk (HT0 z̃0 + HT1 z̃1 ),
(VII.5)
k χk Hχk
where H = HT0 H0 + HT1 H1 . The second crucial observation is that the unique solution vk,r of the local
least squares problem (VII.3) in the (2r)-neighborhood of the vertex k provides a local approximation

15

to the solution x̃ of the global least squares problem (VII.2) in a r-neighborhood of the vertex k ∈ V .
Therefore we can patch vk,r , k ∈ V , together
X
−1 X
r
vr =
χk0
χrk vk,r
(VII.6)
k0 ∈V

k∈V

to generate an approximation to the solution x̃ of the global least squares problem (VII.2) in `p norm,
i.e., there exists δr,σ ∈ (0, 1) such that
kvr − x̃kp ≤ δr,σ kx̃kp
when the radius parameter r ≥ 1 is chosen appropriately. Set
X
−1 X
r
2r −1 2r
J=
χk 0
χrk (χ2r
k Hχk ) χk .
k0 ∈V

(VII.7)

(VII.8)

k∈V

Based on (VII.5), (VII.6) and (VII.7), we propose the following iterative distributed algorithm with initials
z̃0 , z̃1 ∈ `p :

(m−1)
(m−1)

v(m) = J(HT0 z̃0
+ HT1 z̃1
)


 (m)
(m−1)
(m)
z̃0 = z̃0
− H0 v
(VII.9)
(m)
(m−1)

z̃1 = z̃1
− H1 v(m)


 x(m) = x(m−1) + v(m)
for m ≥ 1, where

(0)

(0)

x(0) = 0, z̃0 = z̃0 , z̃1 = z̃1 .

(VII.10)

Remark VII.1. Decompose H = D + R into a diagonal component D and the remainder R. Then the
classical Jacobi method to solve the linear system Hx = HT0 z̃0 + HT1 z̃1 is
x(m) = D−1 (HT0 z̃0 + HT1 z̃1 − Rx(m−1) ), m ≥ 1.

(VII.11)

The above iterative method converges when H is diagonally dominated, which is not necessarily true for
the case in our setting. We observe that for r = 0, the matrix J in (VII.8) is equal to D−1 . Hence the
sequence x(m) , m ≥ 0, in the proposed algorithm (VII.9) and (VII.10) with r = 0 is the same as the one
in the Jacobi method (VII.11) with initial x(0) = 0.
Write Hl = (hl (i, j))i,j∈V and z̃l = (z̃l (i))i∈V , l = 0, 1. For the distributed implementation of the iterative algorithm (VII.9) and (VII.10), each agent k ∈ V is required to transmit information to its neighboring
vertices in B(k, 2r + 2σ), and to store the number mk = µ(B(k, r)) of its neighboring vertices in B(k, r)
e l,k = (hl (i, j))i∈B(k,2r+σ),j∈B(k,2r+2σ) , l = 0, 1.
and four matrices Hl,k = (hl (i, j))i∈B(k,2r+σ),j∈B(k,2r) and H
Shown in Algorithm VII.1 is a distributed implementation of the iterative algorithm (VII.9) and (VII.10),
where every vertex k ∈ V is required to store data of size O((r + σ)2d ), to perform O((r + σ)2d ) algebraic
manipulations in each iteration, and to transmit data to its (2r + 2σ)-neighborhood twice in each iteration.
In the next theorem, we further show that the iterative algorithm (VII.9) and (VII.10) converges
exponentially when r is appropriately selected.
Theorem VII.2. Let 1 ≤ p ≤ ∞, G be a graph satisfying Assumptions II.1 and II.2, (H0 , H1 ) be a graph
filter bank satisfying Assumptions IV.1, IV.2 and IV.4, κ > 1 be as (IV.17), the condition number of the
matrix H := HT0 H0 + HT1 H1 , and let (G0,L , G1,L ) be as in (VI.4). Set

θ 
(D1 (G))2 (2σ + 1)d κ2
exp −
r (3r + 2σ + 1)d ,
(VII.12)
δr,σ :=
κ−1
2σ

16

Algorithm VII.1 Iterative Distributed Reconstruction Algorithm
Inputs: stop criterion ε and observations z̃l,k = (z̃l (i))i∈B(k,2r+σ) for l = 0, 1.
Operation: Compute Fk = HT0,k H0,k + HT1,k H1,k , find its inverse (Fk )−1 , and then compute Gl,L;k :=
(Fk )−1 HTl,k , l = 0, 1.
(0)
(0)
(0)
Initialization: xk = 0, z̃0,k = z̃0,k and z̃1,k = z̃1,k .
Iteration:
(m)
(m)
1) uk = G0,L;k z̃0,k + G1,L;k z̃1,k and write uk = (uk (i))i∈B(k,2r) .
2) Communicate to all P
vertices i ∈ B(k, r)\{k} to send data uk (i) and receive data ui (k).
3) Produce v(k) = m1k i∈B(k,r) ui (k).
4) Communicate to all vertices i ∈ B(k, 2r + 2σ)\{k} to send data v(k) and receive data v(i), and
then generate a vector vk = (v(i))i∈B(k,2r+2σ) .
(m+1)
(m)
(m+1)
(m)
e l,k vk , l = 0, 1.
5) Update xk
= xk + vk and z̃l,k
= z̃l,k − H
(m+1)
6) Evaluate kvk k∞ ≤ ε. If yes, terminate the iteration and output xk
. Otherwise, set m = m + 1.
(m+1)
Outputs: xk
.

where θ = ln(κ/(κ − 1)), σ ≥ 1 is the bandwidth of the analysis filter bank (H0 , H1 ), and d and D1 (G)
are the Beurling dimension and density of the graph G respectively. Take z̃0 , z̃1 ∈ `p , and let x(m) , m ≥ 0,
be as in (VII.9) and (VII.10). If the radius parameter r is so chosen that
δr,σ ∈ (0, 1),

(VII.13)

then x(m) , m ≥ 0, converges to the least squares solution x̃ in (VII.1) exponentially,
kx(m) − x̃kp ≤ (δr,σ )m kx̃kp , m ≥ 0.

(VII.14)

Remark VII.3. For l = 0, 1, we can apply (VII.9) to prove by induction on m that
(m)

z̃l

− (z̃l − Hl x̃) = −H0 (x(m) − x̃), m ≥ 0.
(m)

This together with Theorem VII.2 implies that z̃l , m ≥ 1, in the iterative algorithm (VII.9) and (VII.10)
converges to z̃l − Hl x̃ exponentially, where l = 0, 1.
By (VII.12) and (VII.14) in Theorem VII.2, the iterative distributed algorithm (VII.9) and (VII.10) has
fast convergence rate when a large radius parameter r is chosen. In that case, heavier burden arises at each
iteration, which implies that each vertex in the graph G should have more data storages, better computing
abilities and stronger communication capacities in real world applications. Shown in Tables I and II are
the average Em,r of kx(m) − xk∞ /kxk∞ over 50 trials versus the number m ≥ 1 of iterations and the
spln
radius parameter r ≥ 0, where (Hspln
0,n , H1,n ), n = 1, 2 are used as analysis filter banks, the signal x in
Tables I and II is randomly selected on the Minnesota traffic graph and on the random geometric graph
RGG4096 in Figure 2 respectively. This demonstrates that the iterative distributed algorithm (VII.9) and
(VII.10) converges faster for larger radius r, and the original signal can be well approximated in one step
when a large radius r is chosen, see Tables I and II.
By (VII.12) and Theorem VII.2, there is a radius parameter r0 such that the iterative distributed algorithm
(VII.9) and (VII.10) converges exponentially whenever r ≥ r0 . We can select the above radius parameter
r0 to be independent of the size of the graph G. Our simulation indicates that the iterative distributed
algorithm (VII.9) and (VII.10) with r = 0, i.e. the Jacobi iterative method by Remark VII.1, diverges for
some bounded inputs on the Minnesota traffic graph and on some random geometric graphs, see the first
column of Tables I and II.

17

TABLE I
P ERFORMANCE OF THE ITERATIVE DISTRIBUTED RECONSTRUCTION ALGORITHM TO RECOVER SIGNALS ON THE M INNESOTA TRAFFIC
GRAPH

n=1
Em,r r
m
1
2
3
4
5
10
Em,r r
m
1
2
3
4
5
7
10
14

0
.4155
.1355
.0547
.0226
.0098
.0002

1
.2220
.0238
.0039
.0006
.0000
.0000
n=2

2

3

4

6

.0375
.0007
.0000
.0000
.0000
.0000

.0160
.0001
.0000
.0000
.0000
.0000

.0033
.0000
.0000
.0000
.0000
.0000

.0003
.0000
.0000
.0000
.0000
.0000

0

1

2

3

4

6

1.1988
1.1662
1.4550
1.4697
2.5386
4.7083
12.7921
52.4168

.6563
.2315
.1162
.0567
.0296
.0082
.0013
.0001

.3187
.0518
.0125
.0026
.0006
.0000
.0000
.0000

.1523
.0136
.0017
.0002
.0000
.0000
.0000
.0000

.0725
.0029
.0002
.0000
.0000
.0000
.0000
.0000

.0178
.0002
.0000
.0000
.0000
.0000
.0000
.0000

TABLE II
P ERFORMANCE OF THE ITERATIVE DISTRIBUTED RECONSTRUCTION ALGORITHM TO RECOVER SIGNALS ON THE RANDOM GEOMETRIC
GRAPH RGG4096 IN F IGURE 2
n=1
Em,r r
m
1
2
3
4
10
19
Em,r r
m
1
2
3
4
5
8

0
.4182
.1963
.1143
.0699
.0050
.0001

1
.1301
.0083
.0008
.0000
.0000
.0000
n=2

2

3

4

5

.0156
.0002
.0000
.0000
.0000
.0000

.0035
.0000
.0000
.0000
.0000
.0000

.0006
.0000
.0000
.0000
.0000
.0000

.0001
.0000
.0000
.0000
.0000
.0000

0

1

2

3

4

5

1.5267
2.8586
6.6794
16.4089
40.9430
672.8632

.4674
.1098
.0374
.0121
.0041
.0002

.1487
.0120
.0014
.0002
.0000
.0000

.0437
.0011
.0000
.0000
.0000
.0000

.0159
.0001
.0000
.0000
.0000
.0000

.0049
.0000
.0000
.0000
.0000
.0000

VIII. D ISTRIBUTED DENOISING
Given an NSGFB with analysis filter bank (H0 , H1 ) and synthesis filter bank (G0 , G1 ), we propose a
denoising technique with hard thresholding operator Tτ applied to the high-pass subband signal, where
Tτ (t) = sgn(t)(|t| − τ )+ is the hard thresholding function with threshold value τ ≥ 0, cf. [26], [27], [45],
[46]. Presented in Figure 5 is the block diagram of the proposed denoising procedure. In this section, we
demonstrate the performance of the proposal denoising procedure associated with spline NSGFBs, which
can be implemented in a distributed manner.

18

Fig. 5. Block diagram of the proposed denoising procedure, where x is the noisy input and x̃ is the denoised output.

TABLE III
D ENOISING PERFORMANCE ON THE M INNESOTA TRAFFIC GRAPH MEASURED WITH THE STANDARD `2 -SNR
η
Input `2 -SNR
graphBior
OSGFB
PRT
NSGFB-B1
NSGFB-B2
NSGFB-L1
NSGFB-L2

1/32
34.89
34.43
38.25
35.31
37.50
37.25
38.49
37.25

1/16
28.85
28.91
32.59
29.41
31.45
30.74
32.44
30.67

1/8
22.83
24.06
24.44
23.74
25.43
24.95
26.42
24.91

1/4
16.82
18.21
16.70
18.46
18.91
18.53
19.25
18.16

1/2
10.81
12.79
12.54
15.45
13.18
13.32
13.82
13.33

1
4.75
7.39
4.69
12.77
7.39
7.68
8.34
7.88

In the simulations, the noisy input is
x = xo +  ,

(VIII.1)

where xo = (xo,i )i∈V is the original graph signal and the noise  = (i )i∈V has value i at vertex i ∈ V
spln
randomly selected in the range [−η, η]. The spline NSGFBs have analysis spline filter banks (Hspln
0,n , H1,n )
L,spln
in (III.9) and synthesis spline filter banks being either (GB,spln
, GB,spln
) in (V.8) or (G0,n
, GL,spln
0,n
1,n
1,n ) in
(VI.9), where n ≥ 1. They are abbreviated by NSGFB-Bn and NSGFB-Ln respectively. The denoising
procedure is performed by retaining the low-pass subband signal z0 = Hspln
0,n x and applying the hard threshspln
olding operation Tτ to the high-pass subband signal z1 = H1,n x, where τ > 0 is chosen appropriately.
Thus the denoised output is
x̃ = GB,spln
z0 + GB,spln
Tτ (z1 )
0,n
1,n
for NSGFB-Bn, and
L,spln
x̃ = GL,spln
0,n z0 + G1,n Tτ (z1 )

for NSGFB-Ln respectively, where n ≥ 1. For the above denoising procedure, we use 20 log10 kxo kp /kx−
xo kp to measure the input `p -signal-to-noise ratio (`p -SNR) in dB, and 20 log10 kxo kp /kx̃−xo kp to measure
the output `p -SNR in dB, where 1 ≤ p ≤ ∞.
The Minnesota traffic graph is a test bed for various techniques in signal processing on graphs of
medium size ([8], [20], [22], [26]). The denoising performance of the proposed spline NSGFBs on the
Minnesota graph is presented in Table III, where the original signal xo is the blockwise constant function
in Figure 3, the threshold value τ is selected to be 3η, and the input and output `2 -SNRs are the average
values over 50 trials. Shown also in Table III are the performance comparison with the biorthogonal
graph filter bank (graphBior) in [21], the M -channel oversampled graph filter bank (OSGFB) in [22], and
the pyramid transform (PRT) in [26], where the corresponding output `2 -SNRs are calculated from the
accompanying codes in these references. It indicates that the spline NSGFBs and the OSGFB outperform
other two methods in the small noise scenario, the spline NSGFBs have the best performance in the
moderate noise environment, and the PRT stands out from the rest in the strong noisy case.
Presented in Tables IV and V are the denoising performance of spline NSGFBs and the performance
comparision with the graphBior in [21], the OSGFB in [22], and the PRT in [26] on the random geometric
graph RGGN , where N = 4096, the original signal xo is the blockwise polynomial in Figure 3, the
threshold value τ is selected to be 3η, and the input and output `2 -SNRs in Table IV and the input and

19

TABLE IV
D ENOISING PERFORMANCE ON THE RANDOM GEOMETRIC GRAPH RGG4096
η
Input `2 -SNR
graphBior
OSGFB
PRT
NSGFB-B1
NSGFB-B2
NSGFB-L1
NSGFB-L2

1/32
35.06
33.82
31.69
32.89
37.43
36.65
38.86
36.08

1/16
29.04
28.61
26.37
27.51
31.40
30.63
32.87
29.97

1/8
23.02
23.27
20.79
22.44
25.34
24.89
26.61
24.27

1/4
17.01
18.20
16.40
17.70
19.31
19.37
20.45
19.16

MEASURED WITH THE STANDARD

1/2
10.97
13.21
13.40
14.21
13.47
13.80
14.91
13.92

TABLE V
D ENOISING PERFORMANCE ON THE RANDOM GEOMETRIC GRAPH RGG4096
η
Input `∞ -SNR
graphBior
OSGFB
PRT
NSGFB-B1
NSGFB-B2
NSGFB-L1
NSGFB-L2

1/32
34.90
23.45
20.59
23.24
31.84
26.86
29.28
24.66

1/16
28.88
17.28
14.32
17.25
25.16
20.34
22.70
17.91

1/8
22.85
11.12
6.83
11.27
18.71
14.67
16.19
12.27

1/4
16.83
5.34
0.99
5.43
11.05
8.32
9.28
6.72

1
4.95
8.34
11.13
11.81
7.62
8.25
9.40
9.05

MEASURED WITH THE

1/2
10.81
0.32
−1.99
0.39
6.55
3.55
4.08
0.48

`2 -SNR

`∞ -SNR

1
4.79
−4.00
−2.71
−2.15
2.60
1.67
0.35
−0.52

output `∞ -SNRs in Table V are the average values over 50 trials. It is observed that the spline NSGFBs
proposed in this paper outperform the graphBior, OSBFB and PRT in small and moderate noise scenario,
and that the spline NSGFBs have comparable performance with the rest in the strong noisy case. Also
from Tables IV and V, we see that the differences between the input and output `p -SNRs for p = 2, ∞
are in some range. This confirms the conclusions in Proposition V.2 and Corollary VI.4 that the output
noise is dominated by a multiple of the input noise.
Shown in Figure 6 is the input noise  with η = 1/16 and differences between the original signal xo and
the denoised signal x̃ via the graphBior, OSGFB, PRT and spline NSGFBs, where a random geometric
graph RGG4096 , original signal xo and noise  are the same as in Tables IV and V. It indicates that
all denoising techniques have satisfactory performance inside the same strip where the signal has small
variation, and that the spline NSGFBs proposed in this paper achieve better performance visually on noise
suppression than the other three methods do near the boundary of two adjacency strips where the signal
has large variation.
The proposed NSGFBs can be implemented in a distributed manner and they are beneficial to (local)
noise suppression on graphs of very large scale. Our simulations indicate that for random geometric graphs
RGGN with large size N and 1 ≤ p ≤ ∞, the output `p -SNRs of spline NSGFBs have invisible change
for the same input noise level when the graph size N increases.
A PPENDIX
A. Proof of Proposition III.3
The first inequality follows from (III.7). Now we prove the second inequality. Write A = (a(i, j))i,j∈V ,
and define its Schur norm by


X
X
|a(i, j)|, sup
|a(i, j)| .
kAkS = max sup
i∈V

j∈V

j∈V

i∈V

20

Fig. 6. Plotted on the top from left to right are the input noise  , differences between the original signal xo and the denoised signal x̃
obtained by the graphBior, OSGFB and PRT. Shown at the bottom from left to right are differences between the original signal xo and the
denoised signal x̃ obtained by NSGFB-B1, NSGFB-B2, NSGFB-L1 and NSGFB-L2.

It is well known that the filter bound kAkBp , 1 ≤ p ≤ ∞, of a graph filter A is dominated by its Schur
norm,
(A.1)
kAkBp ≤ kAkS for all 1 ≤ p ≤ ∞.
Then it suffices to prove
kAkS ≤ D1 (G)(σ + 1)d kAk∞ .

(A.2)

For any i ∈ V , we obtain
X

|a(i, j)| =

j∈V

X

|a(i, j)| ≤ kAk∞

ρ(i,j)≤σ

X

1

ρ(i,j)≤σ

≤ D1 (G)(σ + 1)d kAk∞ ,
where the second inequality follows from (II.2). Similarly for any j ∈ V , we have
X
|a(i, j)| ≤ D1 (G)(σ + 1)d kAk∞ .

(A.3)

(A.4)

i∈V

Combining (A.3) and (A.4) completes the proof.
B. Proof of Theorem IV.6
The upper bound estimate (IV.16) follows directly from Proposition III.3 and the observation that
1/2

kH0 kB2 + kH1 kB2 ≤ 2kHkB2 .

(A.5)

Now we prove the lower bound estimate (IV.15). Set
B=I−

H
.
kHkB2

(A.6)

Then B has bandwidth 2σ,
kBkB2 ≤ (κ − 1)/κ,
and
H

−1

= (kHkB2 )

−1

∞
X
n=0

Bn .

(A.7)

(A.8)

21

Write H−1 = (g(i, j))i,j∈V . For κ = 1, we have
H−1 = (kHkB2 )−1 I.

(A.9)

Now we consider the case that κ > 1. Set θ = ln(κ/(κ − 1)), and for i, j ∈ V let n0 (i, j) be the minimal
integer such that 2n0 (i, j) ≥ ρ(i, j)/σ. Then
∞
X

|g(i, j)| ≤ (kHkB2 )−1
≤ (kHkB2 )−1

kBn k∞

n=n0 (i,j)
∞
X

kBknB2

n=n0 (i,j)
−1

≤ (kHkB2 ) κ(1 − κ−1 )n0 (i,j)


θ
≤ kH−1 kB2 exp −
ρ(i, j) ,
(A.10)
2σ
where the first inequality follows from (A.8) and the observation that Bn have bandwidth 2nσ, the second
one is true by (III.7), and the third one holds by (III.8) and (A.7).
From (A.9) we immediately get
kH−1 kBp = kH−1 kB2
(A.11)
if κ = 1, and by (A.1) and (A.10), we have
kH−1 kBp ≤ kH−1 kB2 ×
∞
X
sup
≤
≤

X

 θ

exp −
ρ(i, j)
2σ

i∈V n=0
2nσ≤ρ(i,j)<2(n+1)σ
∞

X
−1
kH kB2 sup
e−nθ µ B
i∈V n=0
∞
X
(n
(2σ)d D1 (G)kH−1 kB2
n=0


i, 2(n + 1)σ − 1
+ 1)d (1 − κ−1 )n

 1 (d)
1−t
d
d+1
−1
≤ d!(2σ) D1 (G)κ kH kB2

≤ (2σ)d D1 (G)kH−1 kB2


t=1−κ−1

(A.12)

if κ > 1. Then
kxkp ≤ kH−1 kBp kHT0 kBp kH0 xkp + kHT1 kBp kH1 xkp



≤ d!(2σ)d D1 (G)κd+1 kH−1 kB2

× kHT0 kBp kH0 xkp + kHT1 kBp kH1 xkp
≤ d!2d (σ + 1)2d (D1 (G))2 κd+1 kH−1 kB2

× kH0 kB2 kH0 xkp + kH1 kB2 kH1 xkp
−1/2

≤ d!2d+1 (σ + 1)2d (D1 (G))2 κd+2 kHkB2

× max kH0 xkp , kH1 xkp ,

−1/2

≤ d!2d+1 (σ + 1)2d (D1 (G))2 κd+2 kHkB2
1
× kH0 xkpp + kH1 xkpp p ,

(A.13)

where the second inequality follows from (A.11) and (A.12), the third holds by Proposition III.3, and the
fourth one is true by (A.5) and (IV.17). This proves (IV.15) and completes the proof.

22

C. Proof of Theorem V.1
By (III.3), (V.1), (V.2) and (V.3), we obtain
G0 H0 + G1 H1
sym
sym
sym
= Q0 (Lsym
G )P0 (LG ) + Q1 (LG )Q1 (LG )
Λ)P0 (Λ
Λ) + Q1 (Λ
Λ)P1 (Λ
Λ))U = UT U = I.
= UT (Q0 (Λ
This completes the proof.
D. Proof of Proposition V.2
Set z0 = H0 x and z1 = H1 x. Then
kx̃ − xkp ≤ kG0 (z0 − Ψ0 (z0 ))kp + kG1 (z1 − Ψ1 (z1 ))kp
≤ (kG0 kBp + kG1 kBp )
≤ D1 (G)(σ̃ + 1)d (kG0 k∞ + kG1 k∞ ),

(A.14)

where the first inequality follows from the perfect reconstruction condition (I.1) for the NSGFB constructed
in Theorem V.1, the second one holds by (V.9), and the last estimate is true by Proposition III.3.
E. Proof of Theorem VI.1
By (VI.4) and (A.10), we have
|gl,L (i, j)| ≤ kH−1 kB2 kHl k∞

X
ρ(k,j)≤σ



θ
ρ(i, k)
exp −
2σ

−1

≤ D1 (G)kH kB2 kHl k∞ (σ + 1)d

θ
θ
ρ(i, j) +
, i, j ∈ V,
× exp −
2σ
2
where l = 0, 1. This proves (VI.6).
F. Proof of Theorem VII.2
Set y(m) = x̃ − x(m) and write y(m) = (y (m) (i))i∈V , m ≥ 0. We claim that
(m)

y(m) = H−1 (HT0 z̃0

(m)

+ HT1 z̃1 ), m ≥ 0.

(A.15)

The above claim holds for m = 0, since
(0)

(0) 

y(0) = x̃ = H−1 (HT0 z̃0 + HT1 z̃1 ) = H−1 HT0 z̃0 + HT1 z̃1
by (VII.1) and (VII.10). Inductively for m ≥ 1, we have
y(m) = y(m−1) − v(m)
(m−1)
(m−1) 
= H−1 HT0 z̃0
+ HT1 z̃1
− v(m)

(m)
(m)
= H−1 HT0 z̃0 + HT1 z̃1

where the first and third equalities follow from (VII.9) and the second equality holds by the inductive
hypothesis. This completes the proof of Claim A.15.
2r −1
Write (χ2r
= (gk (i, j))i,j∈B(k,2r) and
k Hχk )
2r+2σ
2r −1 2r
χrk (χ2r
− χ2r
k ) = (g̃k (i, j))i,j∈V , k ∈ V.
k Hχk ) χk H(χk

(A.16)

23

Following the argument used to prove (A.10), we have


θ
−1
ρ(i, j)
|gk (i, j)| ≤ kH kB2 exp −
2σ
for all i, j ∈ B(k, 2r). By (III.7), (A.16) and (A.17), we obtain
g̃k (i, j) = 0
where either i 6∈ B(k, r) or j 6∈ B(k, 2r + 2σ)\B(k, 2r), and


X
θ
|g̃k (i, j)| ≤ kH−1 kB2 kHk∞
exp −
ρ(i, l)
2σ
l∈B(j,2σ)


θ
≤ D1 (G)(2σ + 1)d κ exp −
r+θ
2σ
where i ∈ B(k, r) and j ∈ B(k, 2r + 2σ)\B(k, 2r).
(m)
(m)
Write vk = (vk (i))i∈V , m ≥ 1, k ∈ V . By (VII.9), (A.15), we have
(m)

χrk (vk

(A.17)

(A.18)

(A.19)

2r −1 2r
− y(m−1) ) = χrk (χ2r
k Hχk ) χk
(m−1)
×H(χ2r+2σ
− χ2r
.
k )y
k

Combining the above equation with (A.18) and (A.19), we get
X
(m)
|vk (i) − y (m−1) (i)| =
g̃k (i, j)y (m−1) (j)
j∈B(k,2r+2σ)


θ
r+θ
≤ D1 (G)(2σ + 1) κ exp −

2σ
X
×
|y (m−1) (j)| , i ∈ B(k, r).
d



(A.20)

j∈B(i,3r+2σ)

This together with (VII.9) implies that
|y (m) (i)| = |v (m) (i) − y (m−1) (i)|
X
1
(m)
≤
|vk (i) − y (m−1) (i)|
µ(B(i, r))
k∈B(i,r)


θ
≤ D1 (G)(2σ + 1)d κ exp −
r+θ

2σ
X
(m−1)
×
|y
(j)|

(A.21)

j∈B(i,3r+2σ)

for all i ∈ V and m ≥ 1. Using the above componentwise estimate, we obtain
ky(m+1) kp ≤ δr,σ ky(m) kp , m ≥ 0.
Iteratively applying the above estimate proves (VII.14).

(A.22)

24

R EFERENCES
[1] I. F. Akyildiz, W. Su, Y. Sankarasubramaniam, and E. Cayirci, “Wireless sensor networks: a survey,” Comput. Netw., vol. 38, pp.
393-422, Mar. 2002.
[2] C. Chong and S. Kumar, “Sensor networks: evolution, opportunities, and challenges,” Proc. IEEE, vol. 91, pp. 1247-1256, Aug. 2003.
[3] J. Yick, B. Mukherjee, and D. Ghosal, “Wireless sensor network survey,” Comput. Netw., vol. 52, pp. 2292-2330, Aug. 2008.
[4] R. Hebner, “The power grid in 2030,” IEEE Spectr., vol. 54, pp. 51-55, Apr. 2017.
[5] C. Cheng, Y. Jiang, and Q. Sun, “Spatially distributed sampling and reconstruction,” Appl. Comput. Harmon. Anal., to be published.
[6] F. Chung and L. Lu, Complex Graphs and Networks (CBMS Regional Conference Series in Mathematics 107), American Mathematical
Society, 2006.
[7] R. Coifman and M. Maggioni, “Diffusion Wavelets,” Appl. Comput. Harmon. Anal., vol. 26, pp. 53-94, July 2006.
[8] 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 Process. Mag., vol. 30, pp. 83-98,
May 2013.
[9] A. Sandryhaila and J. M. F. Moura, “Big data analysis with signal processing on graphs: Representation and processing of massive
data sets with irregular structure,” IEEE Signal Process. Mag., vol. 31, pp. 80-90, Sept. 2014.
[10] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on graphs,” IEEE Trans. Signal Process., vol. 61, pp. 1644-1656, Apr.
2013.
[11] A. Sandryhaila and J. M. F. Moura, “Discrete signal processing on graphs: frequency analysis,” IEEE Trans. Signal Process., vol. 62,
pp. 3042-3054, June 2014.
[12] I. Daubechies, Ten Lectures on Wavelets (CBMS-NSF Regional Conference Series in Applied Mathematics 61), SIAM, 1992.
[13] S. Mallat, A wavelet tour of signal processing: the sparse way, Academic Press, 2009.
[14] M. Vetterli and J. Kovacevic, Wavelets and Subband Coding, Prentice Hall PTR, 1995.
[15] M. Crovella and E. Kolaczyk, “Graph wavelets for spatial traffic analysis,” in Twenty-Second Annual Joint Conf. of the IEEE Computer
and Communications ( INFOCOM 2003), IEEE, vol. 3, San Francisco, USA, 2003, pp. 1848-1857.
[16] M. Gavish, B. Nader, and R. R. Coifman, “Multiscale wavelets on trees, graphs and high dimensional data: Theory and applications
to semi supervised learning,” in Proc. of the 27th Int. Conf. on Machine Learning, Haifa, Israel, 2010, pp. 367-374.
[17] G. Shen and A. Ortega, “Transform-based distributed data gathering,” IEEE Trans. Signal Process., vol. 58, pp. 3802-3815, July 2010.
[18] W. Wang and K. Ramchandran, “Random multiresolution representations for arbitrary sensor network graphs,” in Proc. of 2006 IEEE
Int. Conf. on Acoustics Speech and Signal Processing (ICASSP), IEEE, vol. 4, Toulouse, France, 2006, pp. 161-164.
[19] D. K. Hammod, P. Vandergheynst, and R. Gribonval, “Wavelets on graphs via spectral graph theory,” Appl. Comput. Harmon. Anal.,
vol. 30, pp. 129-150, Mar. 2011.
[20] S. K. Narang and A. Ortega, “Perfect reconstruction two-channel wavelet filter banks for graph structured data,” IEEE Trans. Signal
Process., vol. 60, pp. 2786-2799, June 2012.
[21] S. K. Narang and A. Ortega, “Compact support biorthogonal wavelet filterbanks for arbitrary undirected graphs,” IEEE Trans. Signal
Process., vol. 61, pp. 4673-4685, Oct. 2013.
[22] Y. Tanaka and A. Sakiyama, “M-Channel Oversampled Graph Filter Banks,” IEEE Trans. Signal Process., vol. 62, pp. 3578-3590, July
2014.
[23] V. N. Ekambaram, G. C. Fanti, B. Ayazifar, and K. Ramchandran, “Spline-like wavelet filterbanks for multiresolution analysis of
graph-structured data,” IEEE Trans. Signal Inf. Process. Netw., vol. 1, pp. 268-278, Dec. 2015.
[24] R. M. Karp, “Reducibility among combinatorial problems,” in Complexity of Computer Computations, R.E. Miller, J. W. Thatcher, J.
D. Bohlinger, Ed., 1972, pp. 85-103.
[25] H. Q. Nguyen, and M. N. Do, “Downsampling of signals on graphs via maximum spanning trees,” IEEE Trans. Signal Process., vol.
63, pp. 182-191, Jan. 2015.
[26] D. I. Shuman, M. J. Faraji, and P. Vandergheynst, “A Multiscale Pyramid Transform for Graph Signals,” IEEE Trans. Signal Process.,
vol. 64, pp. 2119-2134, Apr. 2016.
[27] N. Tremblay and P. Borgnat, “Subgraph-based filterbanks for graph signals,” IEEE Trans. Signal Process., vol. 64, pp. 3827-3840,
Aug. 2016.
[28] S. Chen, A. Sandryhaila, and J. Kovačević, “Distributed algorithm for graph signal inpainting,” in 2015 IEEE Int. Conf. on Acoustics,
Speech and Signal Processing (ICASSP), IEEE, Brisbane, Australia, 2015, pp. 3731-3735.
[29] T. T. Doan and L. B. Carolyn, “Distributed primal dual methods for economic dispatch in power networks,” arXiv preprint, arXiv:
1609.06287, Sept. 2016.
[30] D. Baker and A. Ephremides, “The architectural organization of a mobile radio network via a distributed algorithm,” IEEE Trans. on
Commun., vol. 29, pp. 1694-1701, Nov. 1981.
[31] K. Hwang and Z. Xu, “Scalable parallel computers for real-time signal processing,” IEEE Signal Process. Mag., vol. 13, pp. 50-66,
July 1996.
[32] J. Vazquez and C. Malcolm, “Distributed multirobot exploration maintaining a mobile network,” in Proc. of Second IEEE Int. Conf.
on Intelligent Systems, IEEE, Varna, Bulgaria, 2004, pp. 113-118.
[33] Q. Shi, H. Chen, H. Chen, and L. Jiang, “Distributed wireless sensor network localization via sequential greedy optimization algorithm,”
IEEE Trans. Signal Process., vol. 58, pp. 3328-3340, June 2010.
[34] Y. Zhang, T. Cao, S. Li, X. Tian, L. Yuan, H. Jia, and A. V. Vasilakos, “Parallel Processing Systems for Big Data: A Survey,” Proc.
IEEE, vol. 104, pp. 2114-2136, Nov. 2016.
[35] Z. Wang and A. C. Bovik, “Mean squared error: love it or leave it?- A new look at signal fidelity measures,” IEEE Signal Processing
Mag., vol. 98, pp. 98-117, Jan. 2009.
[36] Q. Sun, “Localized nonlinear functional equations and two sampling problems in signal processing,” Adv. Comput. Math., vol. 40, pp.
415-458, Apr. 2014.

25

[37] D. I. Shuman, B. Ricaud, and P. Vandergheynst, “Vertex-frequency analysis on graphs,” Appl. Comput. Harmon. Anal., vol. 40, pp.
260-291, Mar. 2016.
[38] M. Penrose, Random Geometric Graphs (Oxford Studies in Probability 5), Oxford University Press, 2003.
[39] M. S. Kotzagiannidis and P. L. Dragotti, “Sampling and Reconstruction of Sparse Signals on Circulant Graphs - An Introduction to
Graph-FRI,” arXiv preprint, arXiv: 1606.08085, June 2016.
[40] A. Aldroubi, A. Baskakov, and I. Krishtal, “Slanted matrices, Banach frames, and sampling,” J. Funct. Anal., vol. 255, pp. 1667-1691,
Oct. 2008.
[41] C. E. Shin and Q. Sun, “Stability of localized operators,” J. Funct. Anal., vol. 256, pp. 2417-2439, Apr. 2009.
[42] R. Tessera, “Left inverses of matrices with polynomial decay,” J. Funct. Anal., vol. 259, pp. 2793-2813, Dec. 2010.
[43] Q. Sun, “Wiener’s lemma for infinite matrices II,” Constr. Approx., vol. 34, pp. 209-235, Oct. 2011.
[44] C. E. Shin and Q. Sun, “Polynomial control on stability, inversion and powers of matrices on simple graphs,” arXiv preprint, arXiv:
1705.07385, May 2017.
[45] A. Sakiyama and Y. Tanaka. “Oversampled graph Laplacian matrix for graph filter banks,” IEEE Trans. Signal Process., vol. 62, pp.
6425-6437, Dec. 2014.
[46] A. Sakiyama, K. Watanabe, and Y. Tanaka, “Spectral graph wavelets and filter banks With low approximation error,” IEEE Trans.
Signal Inf. Process. Netw., vol. 2, pp. 230-245, Sept. 2016.

