来考虑这样一个线性方程组：$$Ax=b$$其中 $A\in {\mathbb{R}}^{m\times n}，x\in {\mathbb{R}}^n，x\in {\mathbb{R}}^n$。

从图中明显看出这是一个欠定方程组，因为只要 A 的行满秩，该方程组有无穷多组解，否则也可能无解（即 b 不在 A 的列空间内，但我们不考虑这种情况）。

这个方程组是怎么和压缩感知扯上关系的呢？因为 矩阵A 可以看成一个从 n 维空间到 m 维空间的线性映射，显然 $n\gg m$，这是一个压缩映射。

现实的场景是： 采集到的信号 x 是 n 维的，利用压缩变换 A 将原信号压缩成 m 维的 b，由于 $m\ll n$，这将将大大提高信息传播和存储的效率。在这里，我们考虑信号 x 是稀疏的，即 n 个维度中大部分元素为零，只有少量的非零元。

上面这个方程组 $Ax=b$ 的目的就是利用压缩的信号 b，恢复原始信号 x。

如果你认为这就是一个简单的求解线性方程组问题的话，那就大错特错了，因为如前所述，这个方程组有无穷多个解！

实际上，原始信号重建问题对应的是一个约束问题：
$$原始问题：\{\begin{array}{l}\min ||x|{|}_0,\\ s.t.Ax=b\end{array}$$

即 在满足约束 $Ax=b$ 的条件下，经可能地减少 x 中非零元的数目。

不幸的是，上述问题并不是一个凸优化问题，因为 $||\bullet |{|}_0$ 表示非零元个数，不是一个凸函数。

取而代之，我们将优化问题中的目标函数换成 $l_1$ 范数和 $l_2$ 范数来看看，即考虑优化问题：
$$替代问题1：\{\begin{array}{l}\min ||x|{|}_2,\\ s.t.Ax=b\end{array}$$
$$替代问题2：\{\begin{array}{l}\min ||x|{|}_1,\\ s.t.Ax=b\end{array}$$
上述问题可以用现有的凸优化求解器快速求解! 因为 $l_2$范数是凸函数，而替代问题2 可以通过一些变换转换成凸问题。

结果如下：

• 图1 对应原始的稀疏信号；
• 图2 对应在$l_2$范数约束下重建的信号；
• 图2 对应在$l_1$范数约束下重建的信号。

从结果可以看出，$l_2$正则化不能保证稀疏性，而 $l_1$正则化可以！

以下是在 matlab 中调用 CVX 的 mosek 求解器，对上述 $l_2,l_1$ 约束问题求解的代码。

clear all

n = 256;
m = 128;

A = randn(m,n);
u = sprandn(n,1,0.1);
% u = rand(n,1);
b = A*u;

figure(1);
subplot(3,1,1); plot(1:n, u);
title('exact solu');

cvx_solver mosek
cvx_begin 
    variable x(n)
    %minimize( max(norm(x, inf), norm(x,1)/sqrt(n)) )
    %minimize ( max(abs(x)))
    minimize (norm(x))
    subject to
        A*x == b
cvx_end
xl2 = x;

subplot(3,1,2); plot(1:n, xl2);
title('l2 solu');



cvx_begin
    variable x(n)
    minimize( norm(x,1) )
    subject to
        A*x == b
cvx_end
xl1 = x;

subplot(3,1,3); plot(1:n, xl1);
title('l1 solu');

fprintf('\n\nl2 error: %3.2e, l1 error: %3.2e\n', norm(u-xl2), norm(u-xl1));

本文内容参考文再文老师凸优化课程的讲义。

考虑线性系统

寻找解决方案的形式为：

如果数据向量在的列空间中，则上式具有精确解：

如果具有零空间分量，我们将无法满足差分方程，因此我们放宽了要求，并要求尽可能小。

这导致最小二乘问题的定义：

一般最小二乘法解是

有多种解决方案。 例如，正规方程式：

提供解决方案:

正规方程的构造，右边的向量，显然在的列空间中。 实际上，我们给出了组合列向量的方案：

现在，我们可以使用高斯消除和分解之类的工具来解决该问题。

在数值问题中，正规方程可能是一个较差的选择，而和奇异值分解（SVD）等方法较为好。

用SVD下的基础投影，投影在范围空间：

算子将个向量投影到范围空间（列空间，图像）上。

SVD奇异值分解为

这意味着伪逆是

正规方程的解是

MATLAB的代码： 

 function [s, residual] = OLS(A, y, m, err)
% Orthogonal Least Squares [1] for Sparse Signal Reconstruction

% Input
% A = N X d dimensional measurment matrix
% y = N dimensional observation vector
% m = sparsity of the underlying signal

% Output
% s = estimated sparse signal
% r = residual 

% [1] T. Blumensath, M. E. Davies; "On the Difference Between Orthogonal 
% Matching Pursuit and Orthogonal Least Squares", manuscript 2007

if nargin < 4
     err = 1e-5;
end

s = zeros(size(A,2),1);
r(:,1) = y; L = []; Psi = [];
normA=(sum(A.^2,1)).^0.5;
NI = 1:size(A,2);

for i = 2:m+1
    DR = A'*r(:,i-1);
    [v, I] = max(abs(DR(NI))./normA(NI)');
    INI = NI(I);
    L = [L' INI']';
    NI(I)=[];
    Psi = A(:,L);
    x = Psi\y;
    yApprox = Psi*x;
    r(:,i) = y - yApprox;
    if norm(r(:,end)) < err
        break;
    end
end

s(L) = x;
residual = r(:,end);

例子：

% Demo

% Generate Sparse Signal
signal = GenSig(150, 10);
% Create Measurement Matrix
A=matrix_normalizer(randn(100,150));
% Create Measured Signal
y = A*signal;

figure
% reconstruct Signal by OLS
[s1, residual] = OLS(A, y, 10);
stem(signal), hold on, stem(s1,'r+')

结果：

  Justin Romberg  压缩感知：利用凸优化实现稀疏信号的重建的科普

1.利用凸优化算法解决压缩感知问题时，问题通常被分为2类，一类是转化为线性方程组的问题（LPs问题），另一类是转化为二次凸规划的问题（SOCPs问题），其中LPs问题采用路径追踪原始对偶（primal-dual）算法，SOCPs问题采用内点算法中的对数障碍（log-barrier）法

带等式约束的最小化l1范数问题描述如下：

A为M*N的欠定矩阵，x是要采样的信号，b是采样值，大意就是，在方程组条件不足的情况下求解出原始信号x，当x的1范数最小的时候，x差不多就是我们想要的结果了。

L1范数的数学表述形式如下：（本质上来说，就是把矩阵里面的数值求取绝对值然后加起来）

由这个问题引出了著名的基追踪算法（BP算法），算法的原理就是不断寻找l1范数最小的x来解释我们采样到的信号b，如果我们找到了一种足够稀疏的信号（l1范数足够小的信号）

那么我们就认为我们的方程组找到了最合适的解。

2.我们经常能够看到最小化l1范数的误差估计问题，误差估计数学问题描述为如下的问题

A是M*N的满秩矩阵，y是得到的采样值，假设，我们通过信道采样y的时候难免信号会有丢失或者噪声，也就是说我们的得到的y本质上是y=y+e,这里的e就是噪声信号。，同样这个噪声估计的问题也是一个LP问题

3.带有二次约束问题的最小化l1范数问题:

这个问题我们用数学公式表达出来是这样的：

在1问题的基础上我们把原来的等式转为二次约束模型。在这里ε表示一个确定的参数，把原来1问题中的带有等式约束的模型，通过引入松弛变量的方式转为了二次凸规划（SOCP）的问题。这个数学公式的意思是，寻找出一个足够稀疏的信号x满足,在这里e代表一个很小的误差变量

4.有界冗余相关的l1范数最小化：

这种算法也叫Dantzig算法

数学问题描述如下：

我们重新描述了上面的问题，转为上面的公式这里的伽马是一个具体的参数，在1问题里面我们重新引入一种松弛变量，这个数学问题的意思是，找到一种尽可能稀疏的x使得，Ax-b的剩余项与A^*的任何一列最大程度的不相关A^*(Ax-b).

所以问题4中的数学模型是一种LP问题

5.在图像上面引入TV模型，重写上面的公式的表述形式

如果我们要把压缩感知的模型应用在图像上，因为图像的梯度是稀疏的，这里我们用梯度来描述整个图像

定义梯度的操作算子：

分别定义了水平方向和竖直方向的梯度算子这里描述为h和v方向上的梯度

因为整个图像每个点都是梯度的，这里我们把梯度全部加起来（也就全变分）

这就是传说中的（全变分）TV模型，（全变分）TV模型国内研究的人也很多，常用来做去噪和复原的算法，相关论文有很多。

有了上面的定义，我们把上面的数学模型应用在图像复原，上面的每一个数学模型都可以写成SOCP（二次规划）问题

A.带有等式约束的最小化TV模型

B.带有二次规划的最小化TV模型

C.Dantzig TV模型

下一章我将会介绍内点算法。