基于自適應局部搜索的進化多目標稀疏重構方法

劉昊霖 池金龍 鄧清勇 彭 鑫 裴廷睿

1(湘潭大學信息工程學院 湖南湘潭 411105)2(物聯網與信息安全湖南省重點實驗室(湘潭大學) 湖南湘潭 411105)3(湖南理工學院信息科學與工程學院 湖南岳陽 414000)

壓縮感知(compressed sensing, CS)把信號采樣和壓縮相結合，突破了傳統的奈奎斯特采樣定律.壓縮感知理論主要分為3個分支：稀疏表示、觀測矩陣的設計和信號重構[1]，其中信號重構在各領域起著重要作用.信號重構的本質是通過M維觀測信號(M?N)恢復一個長度為N的信號.考慮稀疏信號向量x∈RN，根據壓縮感知理論可知，在理想情況下可以恢復出原始信號x：

(1)

而當觀測過程中存在噪聲n時，觀測向量為y=Ax+n(n∈Rm).眾所周知，直接求解l0范數是一個NP困難問題，可以使用一些方法來求解[2].貪婪算法：包括正交匹配追蹤(orthogonal matching pursuit, OMP)[3]和子空間追蹤(subspace pursuit, SP)[4]、迂回式匹配追蹤(detouring matching pursuit, DMP)[5]等；松弛法：包括基追蹤(basis pur-suit, BP)[6]、基追蹤去噪(basis pursuit denoising, BPDN)[7]、匹配追蹤(matching pursuit, MP)[8]；迭代閾值算法(iterative soft-thresholding algorithm, IST)[9]等.

1) 由于2種局部搜索方法產生2個解，保證解的多樣性，提高了全局解的可信度；

2) 通過比較2個解的優劣性，保證下一個解總是由前最優解產生，提高解的精確度；

3) 由于在進化前期對解進行了批量搜索，獲得了足夠多高質量的解，在后期只需選擇一個優勝的局部搜索方法，既保證了后期解是在最優方法下產生的，又減少了計算成本.

1 相關工作

多目標進化算法主要有非支配排序遺傳算法(non-dominated sorting genetic algorithm, NSGA-II)[14]和基于分解的多目標進化算法(multi-objective evolutionary algorithm based on decomposition, MOEAD)[15]兩大類.NSGA-Ⅱ適合估計連續的帕雷托前沿(Pareto front, PF)，而當PF不連續或非凸時很難找到PF的膝蓋區域(knee region, KR).而且，由于該算法采用擁擠距離來控制種群的多樣性，在尋找KR時會被有些次優解干擾導致定位不準確，最優目標向量在PF上難以呈均勻分布.MOEAD算法將目標函數分解成多個子目標函數，加入規范化技術來處理不同比例的目標函數，從而將偏好整合在PF的局部區域，且可以進化出1組在PF上均勻分布的解.

已有文獻將多目標進化算法應用在稀疏重構問題上，其中文獻[16]提出了多組優化算法來提高在有噪聲環境中的重構精度;文獻[17]提出軟閾值進化多目標算法(soft-thresholding evolutionary multi-objective, StEMO)，通過單次運行產生的1組解來同時優化測量誤差和信號稀疏度,最后在PF[18]的KR上獲得最優解;文獻[19]提出了多目標稀疏模型并且將其用于神經網絡;文獻[20]提出了多目標自步學習模型(multi-objective self-paced learning, MOSPL)用于解決在初始化過程中參數設置的敏感問題;另外，多目標進化算法也在圖像處理領域得到了應用[21-22],文獻[23]提出了多目標稀疏解混模型(multi-objective sparse unmixing, MOSU)用于高光譜圖像解混.

2 多目標問題及求解方法

2.1 基本定義

多目標問題可以描述為

minF(x)=(f1(x),f2(x),…,fm(x))T

(2)

s.t.x∈Ω，

其中，Ω?Rm表示決策空間，x=(x1,x2,…,xn)T為n維決策向量，Ω={x∈Rm|hj(x)≤0,j=1,2,…,m}，hj是連續函數，F:Ω→Rm是由m個目標函數組成，Rm是目標函數空間[24].對于解的支配關系有定義，當且僅當：

(3)

成立時，稱x1支配x2，表示為x1x2.若在解空間中不存在x比x*要好，則稱x*為帕累托最優解，表示為x*x.最優解的目標函數值所組成的集合稱為SPF[25]，表示為

SPF={(f1(x*),f2(x*),…,fm(x*))T|x*∈X*}.

(4)

2.2 切比雪夫分解模型

基于分解的多目標進化方法通常是將前沿面的逼近問題轉化為若干個標量優化問題，一般有3種方分解法，分別是加權法[26]、邊界交集法[27]和切比雪夫法[28].加權法針對凸問題具有更好效果，而邊界交集法必須處理等式約束問題，需要設置懲罰因子，但是懲罰因子的值會影響算法性能.本文采用切比雪夫分解方法，通過調整權重向量來獲得不同的帕累托最優解，其模型可表示為

(5)

3 方法論

3.1 多目標稀疏重構模型

稀疏重構問題可以近似等價為求解式(1)中x的非零項的位置，在傳統方法求解中，通常利用權重參數將測量誤差項和稀疏項聚合成一個函數，但是權重參數的值會影響求解精度.為了解決權重參數對重構精度影響的問題，本文采用多目標進化算法，將測量誤差和信號稀疏度作為2個目標函數:

(6)

3.2 ALSEMO算法描述

算法1.ALSEMO算法.

輸出：解集{x1,x2,…,xN}對應的目標函數組成的集合SEF、目標函數值向量F=(F1,F2,…,FN).

步驟1. 初始化：

步驟1.1. 設SEF=?,k=0.

步驟1.2. 計算出任意2個權重向量之間的歐氏距離，然后獲得任意一個權重向量與之最近的T個權重向量并將其作為鄰居.對任意的i=1,2,…,N，假設B(i)=(i1,i2,…,iT)，λi1,λi2,…,λiT是距離λi最近的T個權重向量.

步驟1.3. 初始化{x1,x2,…,xN}為第1代種群，其目標函數向量Fi=F(xi).

步驟1.4. 初始化參考點集合z*.

步驟2. 循環與更新：

Fori=1,2,…,Ndo

步驟2.2. 執行局部搜索方法：分別將2種局部

步驟2.4. 執行競爭策略選出非支配解：通過比較ya，yb的目標函數值大小，選擇出目標函數值小的解作為非支配解ns.

步驟2.6. 更新相鄰的解和SEF：對任意的j∈B(i)，如果gte(ns|λj,z*)≤gte(xj|λj,z*)，則xj=ns，Fj=F(ns)，移除在SEF中所有被F(ns)支配的向量，如果在SEF中沒有向量可以支配F(ns)，那么就把F(ns)添加到SEF中.

步驟3. 重復步驟2直到k=Maxt；

步驟4. 獲得PF解集{x1,x2,…,xN}；

步驟5. 在PF上用基于角度的方法[18]找到KR并且找出一個最優解.

算法1中初始化種群{x1,x2,…,xN}是由3個步驟完成的：1)隨機產生一個從1～N的整數排列作為x的稀疏支撐集S;2)所有的x都設置為N×1的矩陣;3)產生一個1×K的隨機數矩陣，其元素值在0～1之間，然后將這些元素值分別賦值給支撐集S的1～K個索引，K在本文中表示信號稀疏度.

3.3 差分進化操作

y-r={xi1r+θ×(xi2r-xi3r)， 概率為ω，xi1r， 概率1-ω，

(7)

其中，θ和ω是2個差分進化控制參數.步驟2.2中局部搜索方法(見3.6節)用來提高解的質量.

3.4 修剪策略和高斯變異

(8)

(9)

3.5 競爭策略

在3.2節步驟2.4中，對任意的j=1,2,…,m，如果有fj(ya)

3.6 局部搜索方法

對于稀疏優化問題，為了有效地提高解的質量，在差分進化操作后采用梯度迭代軟閾值(gradient iterative soft-thresholding, GIST)的局部搜索, 這樣有利于促進次優個體向最優方向進化.

GIST是一個范數最小化優化方法，它可以解決問題:

F(x)=f(x)+αg(x)，

(10)

(11)

(12)

3.7 自適應局部搜索方法設計

為了提升解的精度，本文將3.6節中的局部搜索方法結合到進化算法中實現自適應局部搜索方法，一方面可以加速算法收斂速度從而在較短的時間內獲得帕雷托最優解集(Pareto optimal set)，在另一方面，有利于保持解的多樣性，從而在PF有均勻的分布.自適應方法結合了2種局部搜索策略，如圖1所示:

Fig. 1 Adaptive local search method framework圖1 自適應局部搜索方法流程

算法2.基于競爭策略的自適應算法.

輸出：最優解集{x1,x2,…,xN}.

步驟1. 初始化：t1,t2=0,j=1,2,…,m，q=1,2,…,n；

步驟6. 重復步驟3與4直到k=Maxt2且k%2=0.

在差分進化操作后，算法2利用2個局部搜索方法ALS -Ⅰ和ALS -Ⅱ產生2個解，然后比較這2個解的目標函數值選出當前優勝解來更新z*和鄰域信息，最后根據優勝次數選出局部搜索方法產生解.

圖2描述了2種局部搜索算法競爭獲得最優解的流程，其中A={xk,k=0,1,…,n}由ALS -Ⅰ產生，B={sk,k=0,1,…,n}由ALS -Ⅱ產生.它們都滿足式(10)，所以這2個解使得成立：

(13)

Fig. 2 The illustration of the optimal solutiongeneration process圖2 最優解的產生流程示例

根據文獻[19]有xk+1≠xk，sk+1≠sk，所以A和B中解的關系可以描述為

(14)

為了分析xk+1和sk+1的關系，首先需要計算這2個解的目標函數值：

(15)

其中，k=0,1,2,…n.xk+1和sk+1是由2種方法產生的，所以有F(xk+1)≠F(sk+1).除此之外還需要考慮2種情況：

1) 當F(xk+1)

F(xk+1)

(16)

其中，α>0，則有f(xk+1)

2) 當F(sk+1)

如果ALS -Ⅰ在前k次迭代中競爭成功，那么按照算法2，ALS -Ⅱ生成B中的sk+1～sN將被淘汰，不在最優解集中.同理，如果ALS -Ⅱ在前k次迭代中競爭成功，則A中的xk+1～xN將不在最優解集中.具體比較過程如圖2中示例所示，其中假設了一種2組解集可能存在的關系，此時2組解之間的支配關系為

{x1s1,s2x2,x3s3,x4s4,…,sk-1xk-1,xksk,xk+1sk+1,xk+2,xk+3,…,xN-1，xN}，

(17)

則帕雷托最優解集可表示為

{x1,s2,x3,s4,…,sk-1,xk,xk+1,xk+2,xk+3,…,xN-1,xN}，

(18)

其中，sk+1～s1已被淘汰.最后，用基于角度的方法(angle-based method)[18]在帕雷托最優解集上找到最優稀疏解.

4 實驗仿真與分析

實驗分為2部分：1)證明稀疏重構問題存在KR，并且此區域的解能夠使測量誤差和稀疏度達到均衡;2)通過比較9種稀疏重構方法來證明了本文提出的ALSEMO算法具有更高的重構精度.為了與StEMO算法對比，本實驗算法最大迭代次數為200，種群數為100，每個子問題的鄰居個數為20.真實模擬信號xor是隨機產生的稀疏信號，它由3個步驟產生：

1) 隨機選擇非零系數的子集，子集基數代表xor中的非零元素；

2) 非零元素值服從標準正態分布；

3)xor被歸一化為單位長度,測量矩陣A是一個元素值服從標準正態分布的高斯隨機矩陣.

4.1 測試影響KR的因素

4.1.1 測試M對KR的影響

高斯隨機測量矩陣的維度M范圍為600～1 800，間隔為200，N=2 000，x的長度n=2 000，原始信號稀疏度K=100.

Fig. 3 The relations among E, , 圖之間的關系

Fig. 4 The box-plots of KR, and for each M圖4 每個M對應的盒圖

4.1.2 測試δ對KR的影響

Fig. 5 The relations among E, , 圖之間的關系

Fig. 6 The box-plots of KR, and for each δ圖6 每個δ對應的盒圖

4.1.3 測試K對KR的影響

在本組實驗中M=1 200，N=2 000，K的取值范圍為100～300，其間隔為100，δ=0.01.對于每個K，其非零元素的位置是隨機的且元素值服從正態分布.隨著K的增加，稀疏求解問題變得復雜困難，為此，將種群大小和最大迭代次數增加至300，以取得更好的PF和重構效果.圖7顯示了KR的位置隨K變化情況，可以看出KR總是落在K的位置，這意味著KR中解的非零項剛好就是原始信號的稀疏度，KR可以提供一個與真實信號相近的解.對每1組K做10次獨立重復實驗，圖8用盒圖顯示了不同的K對KR變化的統計結果.從圖8中可以看出,KR區域中解的零范數隨著K的變化而變化且其值是近似相等.另一方面，E隨K的增加而增加，這是由于信號越復雜，重構越困難，精度就越低.

Fig. 7 The relations among E, , 圖之間的關系

Fig. 8 The box-plots of KR, and for each K圖8 每個K對應的盒圖

4.2 重構誤差的比較與分析

在本組仿真實驗中，真實信號的產生方法和之前一樣，長度為2 000.將ALSEMO算法和其他8種經典的單目標方法[30]在重構性能上進行了比較，這8種方法分別是：OMP，BP，Homotopy Method (HM)，PFP，L1LS，SpaRSA，FISTA，ALM.對每1組實驗，觀測向量y中加入了服從N(0,0.01)分布的噪聲，圖9～12中每一個數據是經過10次獨立重復實驗后所得的平均重構誤差(average estimation error)，表示為Ea.比較算法中所用的容忍參數為ε=0.3,α=0.02.

Fig. 9 The Ea on different M圖9 不同測量維度M下的重構誤差

Fig. 10 The Ea on different noise level圖10 不同噪聲強度δ下的重構誤差

Fig. 11 The Ea on different sparsity ratio圖11 不同稀疏率KN下的重構誤差

Fig. 12 Performance comparison of ALSEMOwith StEMO圖12 比較ALSEMO和StEMO算法性能

在圖9中，稀疏率KN=0.1，M為600～1 800，間隔為200，其顯示了隨著M的變化Ea的變化情況，可以看出ALSEMO算法的Ea低于其他算法，但是當M在1 000～1 600之間時OMP低于ALSEMO，這是因為容忍參數ε極大地影響了OMP的重構性能，并且這些實驗的選擇值ε=0.3恰好是最佳選擇，但是在ε取其他值時可能表現不佳.在圖10中，KN=0.1，噪聲δ的取值范圍為0.002～0.014, 間隔為0.002.隨著δ的增加，所有算法的Ea都增大，但是ALSEMO算法的Ea始終保持最低，說明ALSEMO具有更好的抗噪性能.圖11顯示了隨著稀疏率KN的增加，所有算法的Ea都在增加，但是ALSEMO算法的Ea始終保持最小.3組實驗證明了本文提出的ALSEMO算法要比其他8種算法的重構精度高，這是由于本文建立了多目標模型，采用多目標進化求解方法避免了傳統求解方法中存在參數敏感的問題，從而提高了重構準確度.

圖12通過比較StEMO算法來證明ALSEMO具有更低的重構誤差.在圖12(a)中，隨著M的增加，ALSEMO和StEMO算法的Ea都在下降，但是ALSEMO總保持較低的Ea，說明ALSEMO具有更高的重構精度.這是由于ALSEMO是以MOEAD為框架，相比于StEMO 它對KR的定位更加準確，并且結合基于l和l1范數的2種局部搜索方法自適應地選擇出最優解，提升了解的精度和收斂速度.在圖12(b)中，兩者的Ea都在增加，但ALSEMO總是低于StEMO，說明ALSEMO具有更好的抗噪性能.在圖12(c)中，Ea隨著稀疏率的增加而增加，但是ALSEMO總低于StEMO，這說明了ALSEMO在信號變得復雜時依然能夠保持較低的重構誤差.

5 總 結