0. 图滤波器
有了上一节中图信号的频率定义之后,接下来我们就可以对图滤波器(Graph Filter)进行定义了。1类比于离散型号处理,在图信号处理中,我们将图滤波器定义为对给定图信号的频谱中各个频率分量的强度进行增强或衰减的操作。假设图滤波器为
H
∈
R
N
×
N
H∈R^{N×N}
H∈RN×N,则
H
H
H:
R
N
→
R
N
R^N→R^N
RN→RN,令输出图信号为
y
\boldsymbol y
y,则:
y
=
H
x
=
∑
k
=
1
N
(
h
(
λ
k
)
x
~
k
)
v
k
\boldsymbol y=H\boldsymbol x=∑_{k=1}^N(h(λ_k)\tilde{x}_k)\boldsymbol v_k
y=Hx=k=1∑N(h(λk)x~k)vk 对比上一节的逆图傅里叶变换的展开向量形式的公式,我们可以清楚地看到增强还是衰减是通过h(λ)项来控制的。为了进一步看清楚H的形式,我们需要对上式进行变换:
y
=
H
x
=
∑
k
=
1
N
(
h
(
λ
k
)
x
~
k
)
v
k
=
[
⋮
⋮
⋯
⋮
v
1
v
2
⋯
v
N
⋮
⋮
⋯
⋮
]
[
h
(
λ
1
)
x
~
1
h
(
λ
2
)
x
~
2
⋮
h
(
λ
N
)
x
~
N
]
=
V
[
h
(
λ
1
)
⋯
0
⋮
⋱
⋮
0
⋯
h
(
λ
N
)
]
[
x
~
1
⋮
x
~
N
]
=
V
[
h
(
λ
1
)
⋯
0
⋮
⋱
⋮
0
⋯
h
(
λ
N
)
]
V
T
x
\begin{aligned}\boldsymbol y&=H\boldsymbol x=∑_{k=1}^N(h(λ_k)\tilde{x}_k)\boldsymbol v_k=\begin{bmatrix}⋮&⋮&⋯&⋮\\\boldsymbol v_1&\boldsymbol v_2&⋯&\boldsymbol v_N\\⋮&⋮&⋯&⋮\end{bmatrix}\begin{bmatrix}h(λ_1)\tilde{x}_1\\h(λ_2)\tilde{x}_2\\⋮\\h(λ_N)\tilde{x}_N\end{bmatrix}\\ &=V\begin{bmatrix}h(λ_1)&⋯&0\\⋮&⋱&⋮\\0&⋯&h(λ_N)\end{bmatrix}\begin{bmatrix}\tilde{x}_1\\⋮\\\tilde{x}_N\end{bmatrix}\\ &=V\begin{bmatrix}h(λ_1)&⋯&0\\⋮&⋱&⋮\\0&⋯&h(λ_N)\end{bmatrix}V^\text{T}\boldsymbol x \end{aligned}
y=Hx=k=1∑N(h(λk)x~k)vk=⎣⎢⎢⎡⋮v1⋮⋮v2⋮⋯⋯⋯⋮vN⋮⎦⎥⎥⎤⎣⎢⎢⎢⎡h(λ1)x~1h(λ2)x~2⋮h(λN)x~N⎦⎥⎥⎥⎤=V⎣⎢⎡h(λ1)⋮0⋯⋱⋯0⋮h(λN)⎦⎥⎤⎣⎢⎡x~1⋮x~N⎦⎥⎤=V⎣⎢⎡h(λ1)⋮0⋯⋱⋯0⋮h(λN)⎦⎥⎤VTx于是得到:
H
=
V
[
h
(
λ
1
)
⋯
0
⋮
⋱
⋮
0
⋯
h
(
λ
N
)
]
V
T
=
V
Λ
h
V
T
H=V\begin{bmatrix}h(λ_1)&⋯&0\\⋮&⋱&⋮\\0&⋯&h(λ_N)\end{bmatrix}V^\text{T}=VΛ_h V^\text{T}
H=V⎣⎢⎡h(λ1)⋮0⋯⋱⋯0⋮h(λN)⎦⎥⎤VT=VΛhVT 相较于拉普拉斯矩阵,
H
H
H仅仅改动了对角矩阵上的值,因此,
H
H
H具有以下形式:
H
i
j
=
0
H_{ij}=0
Hij=0,如果
i
≠
j
i≠j
i=j或者
e
i
j
≠
E
e_{ij}≠E
eij=E。也就是说,
H
H
H矩阵只在对角与边坐标上时才可能取非零值。从算子角度来看,
H
x
H\boldsymbol x
Hx描述了一种作用在每个节点一阶子图上的变换操作。一般来说,我们称满足上述性质的矩阵为
G
G
G的图位移算子(Graph Shift Operator),拉普拉斯矩阵与邻接矩阵都是典型的图位移算子。事实上,本章讲解的所有图信号处理的相关理论,都可以被扩展到图位移算子上,并不局限在拉普拉斯矩阵上。
图滤波器具有以下性质:
(1) 线性:
H
(
x
+
y
)
=
H
x
+
H
y
H(\boldsymbol x+\boldsymbol y)=H\boldsymbol x+H\boldsymbol y
H(x+y)=Hx+Hy;
(2) 滤波操作是顺序无关的:
H
1
(
H
2
x
)
=
H
2
(
H
1
x
)
H_1 (H_2 \boldsymbol x)=H_2 (H_1 \boldsymbol x)
H1(H2x)=H2(H1x);
(3) 如果
h
(
λ
)
≠
0
h(λ)≠0
h(λ)=0,则该滤波操作是可逆的。
我们称
Λ
h
Λ_h
Λh为图滤波器
H
H
H的频率响应矩阵,对应的函数
h
(
λ
)
h(λ)
h(λ)为
H
H
H的频率响应函数,不同的频率响应函数可以实现不同的滤波效果。在信号处理中,常见的滤波器有3类:低通滤波器、高通滤波器、带通滤波器。
图1所示为3类滤波器的频率响应函数:从a图可以看出,低通滤波器只保留信号中的低频部分,更加关注信号中平滑的部分;从b图中可以看出,高通滤波器只保留信号中的高频部分,更加关注信号中快速变化的部分;从c图中可以看出,带通滤波器只保留信号特定频段的成分。
理论上,我们希望实现任意性质的图滤波器,也就是实现任意类型函数曲线的频率响应函数。通过逼近理论,我们可以用泰勒展开——多项式逼近函数去近似任意函数。下面,我们将目光聚焦在拉普拉斯矩阵多项式拓展形式的图滤波器上:
H
=
h
0
L
0
+
h
1
L
1
+
h
2
L
2
+
⋯
+
h
K
L
K
=
∑
k
=
0
K
h
k
L
k
H=h_0 L^0+h_1 L^1+h_2 L^2+⋯+h_K L^K=∑_{k=0}^Kh_k L^k
H=h0L0+h1L1+h2L2+⋯+hKLK=k=0∑KhkLk 其中
K
K
K是图滤波器
H
H
H的阶数。和图信号一样,对于图滤波器,我们也可以同时从空域视角和频域视角来理解。
1. 空域视角
对于
y
=
H
x
=
∑
k
=
0
K
h
k
L
k
x
\boldsymbol y=H\boldsymbol x=∑_{k=0}^Kh_k L^k\boldsymbol x
y=Hx=∑k=0KhkLkx,如果我们设定:
x
(
k
)
=
L
k
x
=
L
x
(
k
−
1
)
\boldsymbol x^{(k)}=L^k\boldsymbol x=L\boldsymbol x^{(k-1)}
x(k)=Lkx=Lx(k−1) 则:
y
=
∑
k
=
0
K
h
k
x
(
k
)
\boldsymbol y=∑_{k=0}^K\boldsymbol h_k\boldsymbol x^{(k)}
y=k=0∑Khkx(k) 通过上式,将输出图信号变成了
(
K
+
1
)
(K+1)
(K+1)组图信号的线性加权。对于式(5),由于
L
L
L是一个图位移算子,因此,
x
(
k
−
1
)
\boldsymbol x^{(k-1)}
x(k−1)到
x
(
k
)
\boldsymbol x^{(k)}
x(k)的变换只需要所有节点的一阶邻居参与计算。总的来看,
x
(
k
)
\boldsymbol x^{(k)}
x(k)的计算只需要所有节点的
k
k
k阶邻居参与,我们称这种性质为图滤波器的局部性。
从空域角度来看,滤波操作具有以下性质:
(1) 具有局部性,每个节点的输出信号值只需要考虑其
K
K
K阶子图;
(2) 可通过
K
K
K步迭代式的矩阵向量乘法来完成滤波操作。
我们以图2为例来看图信号滤波操作在空域的计算过程:
首先将图G的邻接矩阵作为图滤波器:
H
=
A
=
[
0
1
1
0
0
1
0
1
1
0
1
1
0
1
0
0
1
1
0
1
0
0
0
1
0
]
H=A=\begin{bmatrix}0&1&1&0&0\\1&0&1&1&0\\1&1&0&1&0\\0&1&1&0&1\\0&0&0&1&0\end{bmatrix}
H=A=⎣⎢⎢⎢⎢⎡0110010110110100110100010⎦⎥⎥⎥⎥⎤
给定图信号
x
=
[
1
0
0
−
1
−
1
]
T
\boldsymbol x=\begin{bmatrix}1&0&0&-1&-1\end{bmatrix}^\text{T}
x=[100−1−1]T,系数向量
h
=
[
1
0.5
0.5
]
T
\boldsymbol h=\begin{bmatrix}1&0.5&0.5\end{bmatrix}^\text{T}
h=[10.50.5]T,根据公式
y
=
∑
k
=
0
K
h
k
x
(
k
)
\boldsymbol y=∑_{k=0}^K\boldsymbol h_k\boldsymbol x^{(k)}
y=∑k=0Khkx(k),我们可以得到滤波输出的图信号为:
y
=
h
0
x
(
0
)
+
h
1
x
(
1
)
+
h
2
x
(
2
)
\boldsymbol y=h_0 \boldsymbol x^{(0)}+h_1\boldsymbol x^{(1)}+h_2\boldsymbol x^{(2)}
y=h0x(0)+h1x(1)+h2x(2) 其中
x
(
0
)
=
x
=
[
1
0
0
−
1
−
1
]
T
,
x
(
1
)
=
H
x
=
[
0
0
0
−
1
−
1
]
T
,
x
(
2
)
=
H
x
(
1
)
=
[
0
−
1
−
1
−
1
−
1
]
T
。
\begin{aligned}&\boldsymbol x^{(0)}=\boldsymbol x=\begin{bmatrix}1&0&0&-1&-1\end{bmatrix}^\text{T},\\&\boldsymbol x^{(1)}=H\boldsymbol x=\begin{bmatrix}0&0&0&-1&-1\end{bmatrix}^\text{T},\\&\boldsymbol x^{(2)}=H\boldsymbol x^{(1)}=\begin{bmatrix}0&-1&-1&-1&-1\end{bmatrix}^\text{T}。\end{aligned}
x(0)=x=[100−1−1]T,x(1)=Hx=[000−1−1]T,x(2)=Hx(1)=[0−1−1−1−1]T。 带入数据,得到
y
\boldsymbol y
y:
y
=
1
x
+
0.5
x
(
1
)
+
0.5
x
(
2
)
=
[
1
0
0
−
1
−
1
]
+
0.5
[
0
0
0
−
1
−
1
]
+
0.5
[
0
−
1
−
1
−
1
−
1
]
=
[
1
−
0.5
−
0.5
−
2
−
2
]
\begin{aligned}\boldsymbol y&=1\boldsymbol x+0.5\boldsymbol x^{(1)}+0.5\boldsymbol x^{(2)}\\&=\begin{bmatrix}1\\0\\0\\-1\\-1\end{bmatrix}+0.5\begin{bmatrix}0\\0\\0\\-1\\-1\end{bmatrix}+0.5\begin{bmatrix}0\\-1\\-1\\-1\\-1\end{bmatrix}\\&=\begin{bmatrix}1\\-0.5\\-0.5\\-2\\-2\end{bmatrix}\end{aligned}
y=1x+0.5x(1)+0.5x(2)=⎣⎢⎢⎢⎢⎡100−1−1⎦⎥⎥⎥⎥⎤+0.5⎣⎢⎢⎢⎢⎡000−1−1⎦⎥⎥⎥⎥⎤+0.5⎣⎢⎢⎢⎢⎡0−1−1−1−1⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡1−0.5−0.5−2−2⎦⎥⎥⎥⎥⎤
x
(
1
)
\boldsymbol x^{(1)}
x(1)、
x
(
2
)
\boldsymbol x^{(2)}
x(2)的计算在图2中已给出,由于空域滤波操作的局部性,
H
H
H等价于一个聚合邻居的操作算子,如a图、b图中的箭头所示,从a图、b图、c图中可以清楚地看出,
x
(
0
)
→
x
(
1
)
→
x
(
2
)
\boldsymbol x^{(0)}→\boldsymbol x^{(1)}→\boldsymbol x^{(2)}
x(0)→x(1)→x(2)的计算是一个迭代式的过程。
2. 频域角度
H
H
H由于
L
=
V
Λ
V
T
L=VΛV^\text{T}
L=VΛVT,则:
H
=
∑
k
=
0
K
h
k
L
k
=
∑
k
=
0
K
h
k
(
V
Λ
V
T
)
k
=
V
(
∑
k
=
0
K
h
k
Λ
k
)
V
T
=
V
[
∑
k
=
0
K
h
k
λ
1
k
⋯
0
⋮
⋱
⋮
0
⋯
∑
k
=
0
K
h
k
λ
N
k
]
\begin{aligned}H&=∑_{k=0}^Kh_k L^k\\&=∑_{k=0}^Kh_k (VΛV^\text{T})^k\\&=V\bigg(∑_{k=0}^Kh_k Λ^k\bigg) V^\text{T}\\&=V\begin{bmatrix}∑_{k=0}^Kh_k λ_1^k &⋯&0\\⋮&⋱&⋮\\0&⋯&∑_{k=0}^Kh_k λ_N^k \end{bmatrix}\end{aligned}
H=k=0∑KhkLk=k=0∑Khk(VΛVT)k=V(k=0∑KhkΛk)VT=V⎣⎢⎡∑k=0Khkλ1k⋮0⋯⋱⋯0⋮∑k=0KhkλNk⎦⎥⎤ 通过上式,我们可以清楚地看出
H
H
H的频率响应函数为
λ
λ
λ的
K
K
K次代数多项式,如果
K
K
K足够大,我们可以用这种形式去逼近任意一个关于
λ
λ
λ的函数。
如果我们用该滤波器进行类别,则:
y
=
H
x
=
V
(
∑
k
=
0
K
h
k
Λ
k
)
V
T
x
\boldsymbol y=H\boldsymbol x=V\bigg(∑_{k=0}^Kh_k Λ^k\bigg)V^\text{T}\boldsymbol x
y=Hx=V(k=0∑KhkΛk)VTx 上式即为频域视角下的滤波操作,其变换过程由以下3步组成:
(1) 通过傅里叶变换,即
V
T
x
V^\text{T}\boldsymbol x
VTx将图信号变换到频域空间;
(2) 通过
Λ
k
=
∑
k
=
0
K
h
k
Λ
k
Λ^k=∑_{k=0}^Kh_k Λ^k
Λk=∑k=0KhkΛk对频率分量的强度进行调节,得到
y
~
\tilde{\boldsymbol y}
y~;
(3) 通过逆图傅里叶变换,即
V
y
~
V\tilde{\boldsymbol y}
Vy~将
y
~
\tilde{\boldsymbol y}
y~反解成图信号
y
\boldsymbol y
y。
我们假设所有的多项式系数
h
k
h_k
hk构成向量
h
h
h,则
H
H
H的频率响应矩阵为:
Λ
k
=
∑
k
=
0
K
h
k
Λ
k
=
diag
(
Ψ
h
)
Λ^k=∑_{k=0}^Kh_k Λ^k=\text{diag}(\Psi\boldsymbol h)
Λk=k=0∑KhkΛk=diag(Ψh) 其中:
Ψ
=
[
1
λ
1
⋯
λ
1
K
1
λ
2
⋯
λ
2
K
⋮
⋮
⋱
⋮
1
λ
N
⋯
λ
N
K
]
\Psi=\begin{bmatrix}1&λ_1&⋯&λ_1^K\\1&λ_2&⋯&λ_2^K\\⋮&⋮&⋱&⋮\\1&λ_N&⋯&λ_N^K \end{bmatrix}
Ψ=⎣⎢⎢⎢⎡11⋮1λ1λ2⋮λN⋯⋯⋱⋯λ1Kλ2K⋮λNK⎦⎥⎥⎥⎤为范德蒙矩阵,我们可以反解的带多项式系数:
h
=
Ψ
−
1
diag
−
1
(
Λ
k
)
h=\Psi^{-1} \text{diag}^{-1} (Λ^k)
h=Ψ−1diag−1(Λk)其中,
diag
−
1
\text{diag}^{-1}
diag−1表示将对角矩阵编程列向量。
上式说明给定我们想要的频率响应矩阵,可以通过上式反解得到多项式系数。显然,这个十字对于特定性质的图滤波器的设计具有十分重要的意义。
图3是对图信号进行低通滤波的操作示意图。从中可以看到,相较于原始图信号,对其进行低通滤波后,图信号变得更加平滑了。
从频域角度来看,我们可以对图滤波操作作出以下性质总结:
(1) 从频域视角能够更加清晰地完成对图信号的特定滤波操作;
(2) 图滤波器如何设计具有显式的公式指导;
(3) 对矩阵进行特征分解是一个非常耗时的操作,具有
O
(
N
3
)
O(N^3)
O(N3)的时间复杂度,相比空域视角中的矩阵向量乘法而言,有工程上的局限性。
参考文献
[1] 刘忠雨, 李彦霖, 周洋.《深入浅出图神经网络: GNN原理解析》.机械工业出版社. ↩︎