壓縮感知重構算法之基追蹤降噪(Basis Pursuit De-Noising, BPDN)


題目:壓縮感知重構算法之基追蹤降噪(Basis Pursuit De-Noising, BPDN)

        本篇來探討基追蹤降噪(Basis Pursuit De-Noising, BPDN),與基追蹤不同之處在於,基追蹤降噪在模型中考慮到了噪聲的存在,而這在實際中是非常有意義的。因為考慮到了噪聲,所以不同於BP的最優化模型可以轉化為線性規划問題,BPDN的最優化模型可以轉化為二次規划問題。

1、基追蹤降噪的提出

        文獻【1】為了使基追蹤能夠適應有噪聲的數據,提出了基追蹤降噪:

        然后將基追蹤降噪轉化為擾動線性規划(perturbed linear program)問題:

        基追蹤降噪最優化模型中的參數λ可按如下規划選取:

        並沒有查到簡單的MATLAB函數求解擾動線性規划問題。作者發布了文獻【1】中的算法MATLAB源代碼工具箱Atomizer(下載鏈接:http://sparselab.stanford.edu/atomizer/),並有一個說明文檔,在工具箱\Documentation目錄下(AboutAtomizer1208.ps,若想打開需要轉換一下格式,不成功者可以參考鏈接http://download.csdn.net/detail/jbb0523/9583505)

2、利用l1-magic工具箱實現基追蹤降噪

        l1-magic工具箱(鏈接:http://users.ece.gatech.edu/~justin/l1magic/)解決的七類問題中,並沒有出現文獻【1】式(5.1)類型,但在文獻【2】中提到:


        即問題(3)與問題(1)等價。下面是l1-magic的(P2)問題:

        可以看出文獻【2】的問題(3)與l1-magic中的問題(P2)等價,因此,我們完全可以使用l1-magic工具箱中求解(P2)的函數l1qc_logbarrier.m解決基追蹤降噪(使用時注意,函數l1qc_logbarrier還調用了工具箱其它的函數)。

3、將基追蹤降噪轉化為二次規划問題

        文獻【2】給出了將基追蹤降噪轉化為二次規划問題的過程:


        對於線性代數和矩陣分析基礎不好的人,這個過程可能並不是很容易看明白,現詳細推導如下:

Step0:

與基追蹤推導類似,將變量x變為兩個非負變量uv之差:

Step1:

 ,其中 

Step2:

Step3:

因為是一個數(或者說是1×1的矩陣)

所以

Step4:

因為

所以

Step5:

Step6:

綜合Step3、Step4、Step5的結果,Step2可化為:

Step7:

與基追蹤推導類似:

其中, 表示轉置

Step8:

綜合Step6與Step7的結果,Step0可化為:

其中

Step9:

因為是一個與變量無關的常數,因此不影響最優化極值點的位置

所以Step8可寫為標准的二次規划形式:

        求得最優化解z0后可得變量x的最優化解x0=z0(1:n)-z0(n+1:2n)。

4、基於quadprog的基追蹤降噪MATLAB代碼(BPDN_quadprog.m)

        根據以上推導結果,可以寫出基於MATLAB自帶二次規划函數quadprog的代碼:

function [ x ] = BPDN_quadprog( y,A,tao )
%BPDN_quadprog(Basis Pursuit DeNoising with quadprog) Summary of this function goes here
%Version: 1.0 written by jbb0523 @2016-07-22
%Reference:Nowak R D, Wright S J. Gradient projection for sparse reconstruction:
%Application to compressed sensing and other inverse problems[J].
%IEEE Journal of selected topics in signal processing, 2007, 1(4): 586-597.
%(Available at: http://pages.cs.wisc.edu/~swright/papers/FigNW07a.pdf)
[y_rows,y_columns] = size(y);
if y_rows<y_columns
y = y';%y should be a column vector
end
n = size(A,2);
%according to section II-A of the reference
b = A'*y;
c = tao*ones(2*n,1)+[-b;b];
B = [A'*A,-A'*A;-A'*A,A'*A];
lb = zeros(2*n,1);
z0 = quadprog(B,c,[],[],[],[],lb);
x = z0(1:n) - z0(n+1:2*n);
end

5、基追蹤降噪單次重構測試代碼

        測試代碼與OMP測試代碼有改動,首先是測量矩陣Phi使用函數orth進行了正交化,其次是觀測值y加入了噪聲e,可以用比較軟件Beyond Compare具體比較不同之處: 

%壓縮感知重構算法測試  
clear all;close all;clc;
M = 64;%觀測值個數
N = 256;%信號x的長度
K = 10;%信號x的稀疏度
Index_K = randperm(N);
x = zeros(N,1);
x(Index_K(1:K)) = 5*randn(K,1);%x為K稀疏的,且位置是隨機的
Psi = eye(N);%x本身是稀疏的,定義稀疏矩陣為單位陣x=Psi*theta
Phi = randn(M,N);%測量矩陣為高斯矩陣
Phi = orth(Phi')';
A = Phi * Psi;%傳感矩陣
sigma = 0.005;
e = sigma*randn(M,1);
y = Phi * x + e;%得到觀測向量y
%% 恢復重構信號x
tic
lamda = sigma*sqrt(2*log(N));
theta = BPDN_quadprog(y,A,lamda);
x_r = Psi * theta;% x=Psi * theta
toc
%% 繪圖
figure;
plot(x_r,'k.-');%繪出x的恢復信號
hold on;
plot(x,'r');%繪出原信號x
hold off;
legend('Recovery','Original')
fprintf('\n恢復殘差:');
norm(x_r-x)%恢復殘差

運行結果如下:(信號為隨機生成,所以每次結果均不一樣)

         1)圖:

        2)Command Windows

Optimization terminated: relative functionvalue changing by less

 thansqrt(OPTIONS.TolFun), no negative curvature detected in current

 trust region model and the rate of progress(change in f(x)) is slow.

Elapsed time is 8.706162 seconds.

 恢復殘差:

ans =

   0.2723

6、結束語

        從Command Windows輸出可以看到,本次運行花費了8.706162秒,誤差為0.2723。相比於之前的單次重構測試,這個運行時長非常長,誤差非常大。為什么呢?

        有關運行時間長,是不是由於該二次規划變量只有一個非負約束,故可行域范圍比較大,所以時間會很長;而誤差大也是可以理解的,因為本來該方法的觀測值就是含有噪聲的,但問題也就來了,這不是基追蹤“降噪”么?為什么沒把“噪聲”給“降”掉呢?

7、參考文獻:

1Chen S S, Donoho D L, Saunders MA.Atomicdecomposition by basis pursuit[J]. SIAM review, 2001, 43(1): 129-159.(Available at: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.4272&rep=rep1&type=pdf)

 2Nowak R D, Wright SJ. Gradientprojection for sparse reconstruction: Application to compressed sensing andother inverse problems[J]. IEEE Journal of selected topics in signalprocessing, 2007, 1(4): 586-597. (Available at: http://pages.cs.wisc.edu/~swright/papers/FigNW07a.pdf)


注意!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。



 
粤ICP备14056181号  © 2014-2021 ITdaan.com