压缩感知 - 简书

一. 压缩感知理论简介

1.传统采样和压缩感知

1.1 传统采样理论

传统采样理论基于Nyquist采样定理的不足之处

采样速率必须达到信号带宽的两倍以上才能精确重构信号。信息量↑ → 带宽↑ → 硬件成本↑
为了降低存储、处理和传输的成本 ,采用压缩方式以较少的比特数表示信号 ,大量的非重要的数据被抛弃。高速采样再压缩的过程浪费了大量的采样资源。

问题：能否以远低于奈奎斯特采样定理要求的采样速率获取信号 ,同时又可以完全恢复信号 ? 即能否将对信号的采样转变成对信息的采样 ?

1.2 压缩感知理论 $\pmb{^{[1]}}$

压缩感知(Compressed sensing)理论与Nyquist采样定理不同，它指出：

只要信号是可压缩的或在某个变换域是稀疏的 ,那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上 ,然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原信号 ,可以证明这样的投影包含了重构信号的足够信息。

图 3 压缩感知理论框架

其中，待测信号 $x$ 是 $N\times1$ 的列向量，在变换基 $\pmb{\psi}$ 下表示为：
$x=\sum_{k=1}^N\psi_k\alpha_k=\pmb{{\psi}\matrix{\alpha}}\quad\quad\quad\quad\quad\quad\quad (1)$
$x$ 在 $\pmb{\psi}$ 域下是 $K-$ 稀疏的，其中系数 $\alpha_k=<x,\psi_k>$ ，系数向量 $\pmb{\alpha}=\pmb{\psi^T}x$ 。

M次观测下，观测数据 $y$ 是 $M\times1$ 的列向量，与 $\pmb{\psi}$ 不相关的测量矩阵 $\phi$ ： $M\times{N}(M\ll{N})$ 。对系数向量 $\pmb{\alpha}$ 进行线性变换，得到观测集合 $y:M\times1$ 。那么就可以利用优化求解方法从观测集合中精确地重构原始信号。

对于方程 $y=Ax$

传统采样下，为完整恢复待测信号 $x$ ，测量次数 $M$ 至少为 $N$ ,

若 $M<N$ （方程个数小于未知量个数）, $y=Ax$ 方程为欠定方程，无法从 $y$ 的 $M$ 个测量值中解出信号 $x$ 。

若待测信号 $x$ 为 $K-$ 稀疏信号 $K<M\ll{N}$ 时，利用有效恢复算法仍可重建原始信号。

即利用 $0-$ 范数意义下的优化问题求解 $x$ 的精确或近似逼近:
$\min||\pmb{\psi^T}x||_0 \quad\quad\quad{s.t.\quad\quad{y=\phi\psi^Tx}}\quad\quad(2)$
求得信号 $x$ 在 $\pmb{\psi}$ 基上的最稀疏表示 $\alpha$ 。带入式（1）求得 $x$ 。（非凸问题）

式（2）需列出 $x$ 中所有非零项位置的 $C_N^K$ 中可能的线性组合，才能得到最优解，是NP难问题，不可解。 $\pmb{^{[2]}}$

但Chen,Dono、Saunders指出 $\pmb{^{[3]}}$ ，求解一个更加简单的 $l_1$ 优化问题会产生同等的解（ $\pmb{\phi,\psi}$ 不相关），即 $l_1$ 约束也能使得稀疏信号中较小的成分为零，达到稀疏的目的，(2)式等价为：
$\min||\pmb{\psi^T}x||_1 \quad\quad\quad{s.t.\quad\quad{y=\phi\psi^Tx}}\quad\quad(3)$
即就是
$\min||\alpha||_1\quad\quad\quad{s.t.\quad\quad{y=\phi\alpha}}\quad\quad\quad(4)$
满足在测量矩阵 $\phi$ ，使得 $y=\phi\alpha$ 情况下，最小化严格稀疏信号 $\alpha$ 。（凸问题）

2.压缩感知理论主要研究内容

2.1 信号的稀疏表示

稀疏性的定义：

一个实值有限长的 $N\times1$ 维离散信号 $x\in{R^{N\times1}}$ ，可用一组标准正交基 $\pmb{{\psi}}=[\psi_1,\psi_2,..,\psi_k]$ 的线性组合来表示，那么有
$x=\sum_{k=1}^N\psi_k\alpha_k=\pmb{{\psi}\matrix{\alpha}}$
其中 $\alpha_k=<x,\psi_k>$ 构成的 $N\times1$ 维列向量，若 $x$ 在基 $\pmb{\psi}$ 上仅有 $K（K\ll{N}）$ 个非零系数 $\alpha_k$ 时，称 $\pmb{\psi}$ 为信号 $x$ 的稀疏基，称 $x$ 是 $K-$ 稀疏的。 $\pmb{\alpha}$ 是严格 $K-$ 稀疏的。

如何寻找信号的最佳稀疏域，是压缩感知理论的基础和前提。也是信号精确重构的保证。

（1）基函数字典下的稀疏表示：

寻找一个正交基使得信号表示的稀疏系数 $\alpha_k$ 尽可能的少。常用的稀疏变换有，离散傅里叶变换、小波变换、余弦变换等。 $\pmb{^{[4]}}$

（2）冗余字典

冗余函数库来取代基函数，称为冗余字典。从冗余字典中找到具有最佳线性组合的 $K$ 项元素来逼近一个信号，称为信号的稀疏逼近或高度线性逼近。

如何构造合适的冗余字典。

已知冗余字典的前提下，如何使用快速有效方法来稀疏表示信号。

2.2 测量矩阵

测量矩阵 $\phi$ 满足RIP条件：
$(1-\sigma_k)||x||_2^2\le||Ax||_2^2\le(1+\sigma_k)||x||_2^2$
其中， $\sigma_k\in[0,1)$ ,对于任意K-稀疏向量 $x$ 成立。

2.3 信号重构

重构算法四类：

贪婪追踪算法：每次迭代时选择一个局部最优解来逼近原始信号
凸松弛法：将非凸问题转化为凸问题求解找到信号逼近。BP(基追踪)，内点法，梯度投影算法（实验室用），迭代阈值法
组合算法：信号的采样支持通过分组测试快速重建，傅里叶采样,链式追踪和HHS追踪。
深度学习方法

二. GPSR论文

Gradient Projection for Sparse Reconstruction:Application to Compressed Sensing and Other Inverse Problems

1.梯度投影法解决无约束凸问题：

$\min_x\frac{1}{2}\pmb{||y-Ax||}_2^2+\tau||\pmb{x}||_1\quad\quad(5)$

$\pmb{x\in{R}^n,y\in{R}^k},matrix\quad{A}:k\times{n},\tau\ge0$

对于式(3)
$\min||\pmb{\alpha}||_1\quad\quad\quad\quad{s.t.\quad\quad}{y=\phi\alpha}$
$\pmb{\alpha=x}$ 为稀疏信号时，式（3）的解是 $\pmb{x}=0$ ,或者是在某些 $\tau>0$ 情况下，（5）式的解。

关联成像中， $\pmb{x}$ 为像素为 $n$ 原始物体信号， $\pmb{y}$ 为采样次数 $k$ 得到的通探测器信号， $A$ 中每一行为光场数据的行向量化。

2.步骤

将（5）式转化为QP(quadratic program)问题。

令 $\pmb{x=u-v},\pmb{u}\ge0,\pmb{v}\ge0$

$u_i=(x_i)_+,v_i=(-x_i)_+,i=1,2,...,n$

（ $\pmb{u}$ 中元素对应 $\pmb{x}$ 中正数位置，值不变，其他位置为0， $\pmb{v}$ 中元素对应 $\pmb{x}$ 负数位置，值取反，其他位置为0）

$(·)_+$ 定义为： $(x)_+=\max\{0,x\}$

即可将（5）式写为带边界约束的二次规划问题bound-constrained quadratic program (BCQP)
$\min_{u,v}\frac{1}{2}||\pmb{y-A(u-v)}||_2^2+\tau\pmb{1_n^Tu}+\tau\pmb{1_n^Tv}\quad\quad{s.t.\quad\quad{\pmb{u}\ge0,\pmb{v}\ge0}}\quad\quad(6)$
（6）式写为标准形式(BCQP）
$\min_z\quad \pmb{c^Tz}+\frac{1}{2} \pmb{z^TBz}\equiv{F(z)}\quad\quad\quad\quad{s.t.\quad\quad\pmb{z}\ge0}\quad\quad(7)\\ \pmb{z=\begin{bmatrix}u\\v\end{bmatrix},b=A^Ty,c=} \tau\pmb{1}_{2n}+\pmb{\begin{bmatrix}-b\\b\end{bmatrix}，B=\pmb{\begin{bmatrix} A^TA\quad-A^TA\\-A^TA\quad\quad{A^TA} \end{bmatrix}}}$
文中求解（7）式核心：

首先找到一标量参数序列：
$\mathrm{w}^{(k)}=(\mathrm{z}^{(k)}-\alpha^{(k)}\nabla{F(\mathrm{z^{(k)}})})_+$
然后选择第二个标量参数序列：
$\mathrm{z}^{(k+1)}=\mathrm{z}^{(k)}+\lambda^{(k)}(\mathrm{w}^{(k)}-\mathrm{z}^{(k)})$

$\alpha^{(k)}\gt0,\lambda^{(k)}\in[0,1]$

$\mathrm{w}^{(k)}$ 和 $\lambda^{k}$ 的不同产生两种不同的求解算法：
- A. Basic Gradient Projection: The GPSR-Basic Algorithm
  
  每次迭代中沿着负梯度的方向进行搜索，并投影到非负向投影（沿着自变量增。），并执行回溯线搜索Backtracking line search直到函数值获得足够的减少。（BLS算法的思想是，在搜索方向上，先设置一个初始步长，如果步长太大，则缩减步长，直到合适为止。）
  - 初始化沿着负梯度方向能使 $F(z)$ 最小化的步长 $\alpha^{k}$

根据投影定理有：
$-\nabla{F(\mathrm{z^{(k)}})}^T(\mathrm{z^{k+1}-z^k})\ge0$

非负向投影，则 $\mathrm{z^{k+1}-z^k}\ge0$ (正向搜索)，梯度下降方向 $-\nabla{F(\mathrm{z^{(k)}})}^T$ 和搜索方向同号，则 $\nabla{F(\mathrm{z^{(k)}})}\lt0$ 。

防止 $\alpha_0$ 过大或过小，选择 $\alpha_0=mid(\alpha_{min},\alpha_0,\alpha_{max})$ 。（apparently novel）

算法流程

迭代核心代码

for i=1:迭代次数
    while 停止条件
           初始化
           计算导数
        while 1
        % calculate step for this lambda and candidate point，step2的set操作
            du = max(u-lambda*gradu,0.0) - u; 
            u_new = u + du;
            dv = max(v-lambda*gradv,0.0) - v; 
            v_new = v + dv;
            dx = du-dv; 
            x_new = x + dx;

            % evaluate function at the candidate point
            resid_base = A(x_new);
            resid = y - resid_base;
            f_new = 0.5*(resid(:)'*resid(:)) + sum(tau(:).*u_new(:)) + sum(tau(:).*v_new(:));    
            % test sufficient decrease condition
            if f_new <= f + mu * (gradu'*du + gradv'*dv)%step2的判决操作
               %disp('OK')  
            break    %跳出本层while,while keep_going处仍保持运行
            end
        lambda = lambda * lambda_backtrack;        %lambda_backtrack=0.5(前面定义）,放缩为原来的一半
        fprintf(1,'reducing lambda to %6.2e\n', lambda);
      end
      记录结果
      计算停止条件
      ...
    end
end

停止条件

获得充分的高质量的逼近解使困难的，给出能获得较准确的解且避免过度计算

设定 $tolA$
- 目标函数的变化率小于 $tolA$ 时，停止
- $norm(dx(:))/norm(x(:))\lt{tolA}$ 时，停止
- $x$ 中的非零成分数量变化率小于 $tolA$ 时，停止
- 目标函数对于 $u,v$ 求导结果的 $\infty$ 范数变化率小于 $tolA$ 时，停止
B. Barzilai–Borwein Gradient Projection: The GPSR-BB Algorithm

....................................................................................................................

....................................................................................................................

....................................................................................................................