一種改進雅可比算法的頻域臨界采樣圖濾波器組

李威京,蔣俊正

(桂林電子科技大學 信息與通信學院,廣西 桂林 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種改進的雅可比算法近似求解頻域臨界采樣圖濾波器組中的特征矩陣,從而在保證完全重構的前提下達到降低計算復雜度的目的。仿真實驗結果表明,在迭代次數較低時,截斷和并行截斷雅可比算法的去噪性能要優于頂點域采樣圖濾波器組,且接近于傳統頻域臨界采樣圖濾波器組。