一种改进雅可比算法的频域临界采样图滤波器组

李威京,蒋俊正

(桂林电子科技大学 信息与通信学院,广西 桂林 541004)

社交网络、传感器网络及脑神经网络等[1-3]高维不规则数据可建模为图信号。与定义在规则的时间域与空间域的传统离散信号不同,图信号通常定义在不规则的非欧几里得域。传统的信号处理方法不适用于图信号,如何处理这类高维不规则数据已成为一个亟待解决的问题。为此,研究人员提出了图信号处理框架[4],该框架将传统信号处理方法与图论相结合,为处理这类高维不规则数据提供了强有力的工具。其中,由传统滤波器组延伸而来的图滤波器组[5-7]因其稀疏特性和多分辨分析能力而受到了广大研究人员的青睐。

目前,图滤波器组的设计方法主要有顶点域采样图滤波器组和频域采样图滤波器组。文献[8]提出了完全重构的两通道顶点域临界采样图滤波器组,该滤波器组满足正交特性,但不满足紧支撑性。文献[9]提出了双正交图小波滤波器组,该滤波组既满足完全重构特性,又满足紧支撑性。但上述2种顶点域采样滤波器组本质上仅适用于二分图,对于非二分图则需要进行近似处理。文献[10]提出的样条图小波滤波器组满足完全重构特性与图不变性,其对所有拓扑结构的图信号都适用。然而,这几类顶点域采样图滤波器组还存在如下局限性:1)需要选择合适的采样集,才能确保其完全重构;2)其完全重构的采样集不唯一,不同的采样集会影响图滤波器组的整体性能。文献[11]提出了两通道频域临界采样图滤波器组,该频域采样图滤波器组克服了顶点域采样图滤波器的缺点,其完全重构的采样集是唯一的,且对于任意拓扑结构的图信号都满足完全重构特性;但由于其采样操作在频域上进行,需要借助特征分解来获取图模型的特征向量矩阵,这导致了计算复杂度过高。

针对频域临界采样图滤波器组存在的问题,可采用改进的雅可比算法(截断雅可比算法和并行截断雅可比算法)来近似求解拉普拉斯矩阵的特征矩阵。改进的雅可比算法是一种贪婪算法,其每步迭代所求得的稀疏正交矩阵都属于Givens旋转矩阵,将每步迭代求得的Givens旋转矩阵相乘,可得到近似特征矩阵,从而近似求得图信号的频域表示。改进雅可比算法在满足完全重构的条件下降低了频域临界采样图滤波器组的计算复杂度。

1 基础知识

1.1 图信号处理的基础知识

无向图G=(V,E)是由顶点集V={0,1,…,N-1}与边集E组成,其中节点数N=|V|表示图的大小。图G可用权重矩阵W∈RN×N表示,权重矩阵的元素wij表示顶点i与顶点j之间的关联程度,若wij=0,则表示顶点i与顶点j之间没有边相连。根据权重矩阵W,可定义度矩阵D、拉普拉斯矩阵L,度矩阵D是对角矩阵,其对角线元素等于权重矩阵的行和,即,而拉普拉斯矩阵L定义为度矩阵与权重矩阵之差,即L=D-W。对拉普拉斯矩阵L进行特征分解,有

其中,对角矩阵Λ∈RN×N的对角元素由L的特征值组成,其对角元素按从小到大的顺序排列,即0=λ0≤λ1≤…≤λN-1,特征矩阵U=[u0,u1,…,uN-1]的列向量对应于L的特征向量。

图信号可表示为N维向量f∈RN,即f=[f0,f1,…,fN-1]T,其中第i个元素对应于图上第i个节点的信号值。与传统离散信号类似,图信号也有对应的频域表示,图信号f的图傅里叶变换定义为

逆图傅里叶变换定义为

式(2)、(3)称为图傅里叶变换对。

图滤波可从顶点域和频域2个角度考虑。在顶点域,图滤波可表示为拉普拉斯矩阵的多项式,图信号的图滤波过程可表示为

图信号在经过图滤波后,每个顶点的信号值都是其p-hop邻域的信号值的线性组合。在频域,图信号的图滤波过程可表示为

为图滤波器的频率响应,H(λi)为图频率的多项式。

1.2 两通道频域临界采样图滤波器组的结构

图1为两通道频域临界采样图滤波器组的基本结构[11],其中:

图1 两通道频域临界采样图滤波器组

两通道频域临界采样图滤波器组的输入输出关系为

根据输入输出关系,当子带滤波器的频率响应满足

传递函数T=c2I,此时该滤波器组满足完全重构条件[11]。

关于G0(λi)、H0(λi)的设计,现有的设计方法主要有正交[8]和双正交[9]2种类型。两通道频域采样图滤波器组的正交设计取决于H0(λi),而H0(λi)的设计借鉴了经典信号处理中的设计方法[8],即

其余滤波器H1(λi)、G0(λi)、G1(λi)都可由H0(λi)求得,

与此同时,为了确保图滤波器组的完全重构特性,H0(λi)需满足:

对于双正交频域采样图滤波器组而言,其高通滤波器则由低通滤波器定义[9],即

G0(λi)、H0(λi)可通过频谱分解获得。该设计方法类似于经典信号处理中的双正交小波变换[13]。在上述条件的约束下,双正交频域采样图滤波器组的完全重构条件[9]为

2 频域临界采样图滤波器组的快速实现方法

频域临界采样图滤波器组需要计算图信号的图傅里叶变换,因此需对拉普拉斯矩阵L进行特征分解,求得其特征矩阵U,该操作的计算复杂度为O(n3)。当图的规模较小时,可快速求得其特征矩阵,但当图的规模较大时,计算特征矩阵U的代价会变得较大。为降低求解特征矩阵U的计算复杂度,可采用改进的雅可比算法来近似求解特征矩阵U[14]。

2.1 近似求解特征矩阵的问题归结

近似求解特征矩阵U的目的是找到一个近似的特征矩阵,使得以下最优化问题取得最小值:

为了衡量近似特征矩阵与特征矩阵U间的计算复杂度,给出相对复杂度(relative complexity gain,简称RCG)的定义[14]。相对复杂度是特征矩阵U的非零元素个数与近似特征矩阵中的稀疏正交因子Sj的非零元素个数之和的比率,即

其中‖U‖0为矩阵U的0范数。

2.2 截断雅可比算法

截断雅可比算法是雅可比算法[15]的改进,雅可比算法在近似对角矩阵Lj达到某个精度后停止迭代,而截断雅可比算法则预先设置了迭代次数。对于式(12)的最优化问题,若将稀疏正交矩阵Sj约束在Givens旋转矩阵集RG中,并预先设定稀疏正交矩阵Sj的个数J,即可得到截断雅可比算法[14]。n维向量的Givens旋转固定其中n-2个坐标,其余2个坐标则旋转一定的角度,n维Givens旋转矩阵可表示为

其中:p、q为旋转坐标;θ∈[0,2π]为旋转角;c=cosθ;s=sinθ[16]。显然,Givens旋转矩阵由p、q、θ三个参数所决定。截断雅可比算法是一种贪婪算法,其每步迭代都旨在寻找一个Givens旋转矩阵,使得代价函数下降最快[14],即

截断雅可比算法如下:

输入:拉普拉斯矩阵,Givens旋转次数。

输出:近似特征矩阵,近似特征值。

初始化:Lj=L,j=1。

5) 若err_θ1＜err_θ2,则旋转角θ=θ1;若err_θ1＞err_θ2,则旋转角θ=θ2;

6) 计算c=cosθ,s=sinθ及Givens旋转矩阵;

7) 计算Givens旋转后的拉普拉斯矩阵Lj+1=,若j＜J,则返回步骤1),继续迭代;若j=J,则终止迭代,且输出近似特征矩阵=S1S2…SJ,近似特征值=diagLj+1。

2.3 并行截断雅可比算法

式(13)给出的相对复杂度增益只是理论上的度量指标,在实际运算过程中,使用近似特征矩阵运算的时间开销与其相对复杂度增益并不完全成比例。这是因为,在编程环境中(如MATLAB),近似特征矩阵=S1S2…SJ与向量相乘是一个串行过程,而特征矩阵U与向量相乘则是一个并行过程[14]。因此,即使相对复杂度增益较大,使用近似运算的时间也可能相对较长。

由截断雅可比算法改进而来的并行截断雅可比算法[14]较好地解决了上述问题。在每步迭代过程中,截断雅可比算法只做了1次Givens旋转,而并行截断雅可比算法则可做n/2次Givens旋转。因此,对于需要进行J次Givens旋转的近似过程,并行截断雅可比算法只需选择K=2J/n个旋转因子,其中每个旋转因子Sk都是由n/2个不相交Givens旋转所组成的矩阵,其数学表达式可表示为

与截断雅可比算法类似,并行截断雅可比算法也是通过寻找矩阵Lk中绝对值最大的元素来确定Givens旋转,不同的是,所选出的n/2个元素必须保证两两之间不处在同一行或同一列。在每步迭代过程中,虽然并行截断雅可比算法得到的解不是最优解,但其实际运行时间要远低于截断雅可比算法[14]。并行截断雅可比算法如下:

输入:拉普拉斯矩阵L,Gives旋转次数J,每次迭代的Givens旋转次数t=n/2。

输出:近似特征矩阵,近似特征值。

初始化:Lk=L,k=1。

1) 取拉普拉斯矩阵Lk的上三角元素,并按从大到小的顺序进行排列;

2) 取前n/2个最大值元素,且保证两元素不在同一行或同一列;

3) 根据取出的n/2个元素,计算其对应的Givens旋转角,

4) 计算ci=cosθi,si=sinθi,并将其填入矩阵Sk相应的位置;

5) 若err_q1＜err_q2,则旋转角θ=θ1;若err_q1＜err_q2,则旋转角θ=θ2;

6) 计算c=cosθ,s=sinθ及Givens旋转矩阵;

7) 计算Givens旋转后的拉普拉斯矩阵Lj+1=,若j＜J,则返回步骤1),继续迭代;若j=J,则终止迭代,且输出近似特征矩阵=S1S2…SJ,近似特征值=diagLj+1。

2.4 计算复杂度分析

2.4.1 截断雅可比算法

寻找矩阵Lj中绝对值最大的元素是截断雅可比算法中计算复杂度最高的操作,其计算复杂度为O(‖Lj‖0),最坏情况下,该操作的计算复杂度为O(n2)。Givens旋转前后的矩阵Lj与Lj+1仅有p、q两行与p、q两列的元素不相同,其余位置的元素都相同。若当前迭代步骤中所选Givens旋转坐标p、q与前次选择的不同,计算复杂度就能降到O(n)。实际计算过程中,绝大多数情况是2次选择的坐标不相同,因此截断雅可比算法的平均计算复杂度为O(n2+nJ)。求得近似特征矩阵后,该矩阵与向量或矩阵相乘的计算复杂度为O(J)。

2.4.2 并行截断雅可比算法

对矩阵Lk元素进行排序是并行截断雅可比算法中计算复杂度最高的操作。在每步迭代中,最多需要对n(n-1)/2个元素进行排序,该操作的计算复杂度为O(n2logn)。由于并行截断雅可比算法只需迭代K=O(J/n)次,因此整体计算复杂度为O(nJlogn)。在求得近似特征矩阵后,运用该矩阵与向量或矩阵相乘的计算复杂度为O(J)。并行截断雅可比算法中每个旋转因子Sk包含了n/2个Givens旋转,是一个并行计算过程。

为了使上述2种近似对角化算法的计算复杂度低于精确对角化,需要确保迭代次数J＜O(n2),通常迭代次数设为J=O(nlogn)。

3 仿真结果与分析

3.1 图信号去噪

将截断和并行截断雅可比算法与文献[8]、[9]、[11]的算法在图信号去噪上的性能进行对比,所有仿真均在相同环境下进行。用GSPBOX[18]构建了Random sensor、Swiss roll、Community三种图,图上信号由图拉普拉斯矩阵前10%的特征向量求和组成,图2为节点数为128的3种图信号。实验中,图信号所添加的噪声是均值为0、标准差为σ的加性高斯噪声,图滤波器组中的低通信道的阈值设为0,高通信道的阈值设为3σ。

图2 3种不同拓扑结构的图信号

表1～3为不同算法对含噪图信号去噪后的信噪比,去噪结果是50次随机实验结果的平均值。通过对比实验结果,可得如下结论:

表1 sensor图信号在不同算法下去噪后的信噪比 dB

表2 roll图信号在不同算法下去噪后的信噪比 dB

表3 community图信号在不同算法下去噪后的信噪比 dB

1)在Givens旋转次数相同的情况下,截断雅可比算法的去噪效果总体上优于并行截断雅可比算法。这是因为在每步迭代过程中,截断雅可比算法获得的解都是局部最优解,而并行截断雅可比算法获得的解并非局部最优解,这导致了并行截断雅可比算法的收敛速度要慢于截断雅可比算法。

2)随着Givens旋转次数增加,截断雅可比算法与并行截断雅可比算法的去噪效果更接近于文献[11],原因在于,随着迭代次数的增加,2种算法求得的近似特征矩阵不断趋近于精准的特征矩阵。

3)在迭代次数较低时,截断和并行截断雅可比算法取得的去噪效果要优于文献[8-9]。在顶点域采样图滤波器组中,图信号与噪声信号未得到有效分离,2种信号混杂程度较高,从采样图信号中难以有效滤除噪声信号。而对于频域采样图滤波器组而言,由于噪声信息大多分布在高频部分,图信号与噪声信号从频域角度得到了有效分离。因此,频域采样图滤波器组[11]的去噪效果要优于顶点域采样图滤波器组[8-9],而改进雅可比算法是对文献[11]的一种近似,其去噪效果本质上仍优于文献[8-9]。

3.2 图像去噪

在相同噪声环境下对比不同算法对图像的去噪性能。图像可建模为一张图,文献[19]给出了图像的8邻域图表示,其像素对应于图节点,像素值则对应于图信号,每个像素与其上下左右4个像素及2条对角线上的4个像素相连接,其边的权重值由相邻节点的像素值决定[20],如图3所示。

实验选用了如图4(a)所示的硬币图像,并以峰值信噪比(PSNR)和结构相似性(SSIM)作为图像去噪的性能评价指标。图4为不同算法去噪前后的图像,表4、5分别为图像去噪前后的峰值信噪比和结构相似性。从图4和表4、5可看出,与文献[8-9]相比,改进雅可比算法的去噪效果更好;与文献[11]相比,去噪效果稍差,但随着Givens旋转次数的不断增加,改进雅可比算法的去噪效果不断趋近于文献[11]的去噪效果。另外,在低噪声水平下,部分图滤波器组未发挥作用,表现在去噪后信号信噪比低于含噪前,原因在于图滤波器组滤掉的高频部分不仅含有噪声信息,还含有图像的边缘信息,去噪过程中原本的图像信息遭到破坏,而真正的噪声信息未得到有效处理。

表4 图像在不同算法下去噪前后的PSNR dB

图4 硬币图像在不同算法下去噪前后的图像

4 结束语

拉普拉斯矩阵的特征分解是频域临界采样图滤波器组中计算复杂度最高的运算,当图的规模很大时,计算代价过大。针对该问题,使用2种改进的雅可比算法近似求解频域临界采样图滤波器组中的特征矩阵,从而在保证完全重构的前提下达到降低计算复杂度的目的。仿真实验结果表明,在迭代次数较低时,截断和并行截断雅可比算法的去噪性能要优于顶点域采样图滤波器组,且接近于传统频域临界采样图滤波器组。