压缩感知重构算法之Gradient Descent with Sparsification(GraDeS)

来源：互联网发布：网络在线教育平台编辑：程序博客网时间：2024/06/06 12:55

题目：压缩感知重构算法之Gradient Descent with Sparsification(GraDeS)

本篇介绍压缩感知重构算法GraDes(GradientDescent with Sparsification)，该算法是在文献【GargR, Khandekar R. Gradient descent with sparsification: an iterative algorithmfor sparse recovery with restricted isometry property[C]//Proceedings of the26th Annual International Conference on Machine Learning. ACM, 2009: 337-344】中提出的，该论文可由链接【1】下载，另外，该论文在发表的会议上有一个talk，PPT见链接【2】，两位作者的主页参见链接【3】。文章的署名单位为IBM的T.J. Watson Research Center，目前第一作者Rahul Garg已到印度理工学院(德里校区)工作。

1、算法

先说一下符号约定，算法GraDeS是为了求解如下问题：

接下来，我们来看一下算法GraDeS的具体内容：

其中：

函数Ψ(x)定义如下：

符号H_s(·)的含义如下：

也就是说，H_s(x)是将一个n维向量x中最大的s个值保留，其余全部置零，也就是将n维向量x强制变为s稀疏的。文章中紧跟表格Algorithm1有一段有关H_s(·)的描述：

参数γ的取值参见Theorem 2.1：

即若要恢复一个s稀疏的向量x，γ等于1加上矩阵Φ的2s阶有限等距常数δ_2s。另外，在Theorem2.1中要求矩阵Φ的2s阶有限等距常数δ_2s<1/3，还给出了具体的迭代次数公式。但个人感觉，Theorem 2.1的结果并没有实际意义，无论是参数γ还是迭代次数，原因很简单，因为在实际当中矩阵Φ的2s阶有限等距常数δ_2s根本无法得到。

2、GraDeS算法MATLAB代码

作者给出了该函数的MATLAB代码，参见链接【4】，官方版本略微复杂，为了便于理解算法，这里根据Algorithm1，给出一个简化版的GraDeS的实现代码，该算法是如此的简单，以至于我们可以很轻松地写出来：

[plain] view plain copy
 
function [ x ] = GraDeS_Basic( y,Phi,s,gama,epsilon,loopmax )  
%GraDeS_Basic Summary of this function goes here  
%Version: 1.0 written by jbb0523 @2016-07-28  
%Reference:Garg R, Khandekar R. Gradient descent with sparsification: an iterative   
%algorithm for sparse recovery with restricted isometry property[C]//Proceedings of   
%the 26th Annual International Conference on Machine Learning. ACM, 2009: 337-344  
%Available at:http://people.cse.iitd.ernet.in/~rohitk/research/sparse.pdf  
%   Detailed explanation goes here  
    if nargin < 6  
        loopmax = 3000;  
    end  
    if nargin < 5    
        epsilon = 1e-3;    
    end   
    if nargin < 4    
        gama = 1+1/3;    
    end   
    [y_rows,y_columns] = size(y);    
    if y_rows<y_columns    
        y = y';%y should be a column vector    
    end  
    n = size(Phi,2);  
    x = zeros(n,1);%Initialize x=0  
    loop = 0;  
    while(norm(y-Phi*x)>epsilon && loop < loopmax)  
        x = x + Phi'*(y-Phi*x)/gama;%update x  
        %the following two lines of code realize functionality of Hs(.)  
        %1st: permute absolute value of x in descending order  
        [xsorted inds] = sort(abs(x), 'descend');  
        %2nd: set all but s largest coordinates to zeros  
        x(inds(s+1:n)) = 0;  
        loop = loop + 1;  
    end  
end  

注：本来就单纯想按照论文的Algorithm 1中的流程写了个简简单单的函数，但在while循环条件里面仍然加入了最大循环次数的限制，这是因为仿真过程中发现程序经常处于死循环状态。

这里也给出作者的官方版本（参见链接【4】）：

[plain] view plain copy
 
function [x i] = grades(y, phi, s, tol, niters, weight)  
if (nargin < 6)  
    weight = 1 + 1/3;  
end;  
if (nargin < 5)  
    niters = 3e3;  
end;  
if (nargin < 4)  
    tol = 1e-10;  
end;  
if (nargin < 3)   
    disp('Error:  too few arguments');  
    return;  
end;  
  
y2 = sum(y.*y);  
r = y;   
r2 = y2;  
r_last = r2 + 1;  
n = size(phi, 2);  
x = zeros(n, 1);  
i = 0;  
while ((r2 / y2 > tol) && (i < niters) && (r2 / y2 < 1e5) && (r2 < r_last * (1+0.1)))  
    x = x + phi' * r / weight;  
    [xsorted inds] = sort(abs(x), 'descend');  
    x(inds(s+1:n)) = 0;  
    r = y - phi * x;  
    r_last = r2;  
    r2 = sum(r.*r);  
    i = i+1;  
end;  
end  

3、单次重构代码

[plain] view plain copy
 
%压缩感知重构算法测试      
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      
y = Phi * x;%得到观测向量y    
%% 恢复重构信号x      
tic      
theta = GraDeS_Basic(y,A,K);      
% [theta,i] = grades(y,A,K);  
% fprintf('\ni=%d\n',i);  
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：Phi必须经过“Phi= orth(Phi')';”处理，否则根本没有重构结果输出。

4、为什么能够重构？

该算法比以往介绍的任何算法都要简单，但为什么执行以下迭代

最终就可以得到稀疏解x呢？我们来看文章2.1节的一段描述：

也就是说，每次迭代，我们都会根据参数γ沿着目标函数Ψ(x)梯度相反的方向走一步，为了保证每次迭代后仍是s稀疏的，对迭代结果执行硬阈值(hard-thresholding)，也就是函数H_s(·)。

那么这里有两个问题：

一是为什么沿着目标函数Ψ(x)梯度相反的方向走一步会使目标函数变小？

二是执行硬阈值后会不会使目标函数Ψ(x)很大幅度地增加？

对于第一个问题：

实际上，这个类似于最优化方法中的梯度下降法（或称最速下降法），我们举一个简单例子，比如f(x)=x²，求它的梯度（也就是导数）▽f(x)=2x，对于变量x来说只有两个移动方向，要么是正方向，要么是负方向：

假设我们的初始点在x=4，这时梯度▽f(4)=2×4=8>0，梯度方向（单位向量可以通过▽f(x)/|| ▽f(x)||得到）为正方向（即由0指向1），其反方向为负方向（即由0指向-1），从f(x)=x²的图像中可以很容易的知道，从x=4往负方向走一步函数值会下降（当然步子不能迈得太大）；假设我们的初始点在x=-4，这时梯度▽f(4)=2×(-4)=-8<0，梯度方向为负方向（即由0指向-1），其反方向为正方向（即由0指向1），从f(x)=x²的图像中可以很容易的知道，从x=-4往正方向走一步函数值会下降（当然步子也不能迈得太大）。因为这个自变量只有一个，所以可移动的方向只有两个，即在一条直线上移动；如果自变量有两个的话，则移动范围是在一个平面上。求解只有一个变量函数的极值问题即为一维极值问题，求解有多个变量函数的极值问题即为多维极值问题。

对于第二个问题：

为了解释第二个问题，先来看一段文章内容：

注：w.r.t.意思是“关于”，是英文“with regard to”或“withreference to”或“with respect to”的简称。

这段内容首先是对该算法的原理进行了简要的描述，这个类似于分析第一个问题中的引用的那段话；紧跟着有一句非常关键的话“the RIP of Φ impliesthat the sparsification step does not increase the error Ψ(x) by too much”，也就是说Φ的RIP保证了稀疏步骤(sparsificationstep，即执行硬阈值步骤)不会使Ψ(x)误差增加很多。但刨根问底，为什么呢？接下来提到这篇论文重要的贡献是分析了硬阈值函数H_s(·)对Ψ(x)的影响，也就是Lemma 2.4。

引理2.4是由下面内容引出的：