来考虑这样一个线性方程组:
其中
。
从图中明显看出这是一个欠定方程组,因为只要 A 的行满秩,该方程组有无穷多组解,否则也可能无解(即 b 不在 A 的列空间内,但我们不考虑这种情况)。
这个方程组是怎么和压缩感知扯上关系的呢?因为 矩阵A 可以看成一个从 n 维空间到 m 维空间的线性映射,显然 ,这是一个压缩映射。
现实的场景是: 采集到的信号 x 是 n 维的,利用压缩变换 A 将原信号压缩成 m 维的 b,由于 ,这将将大大提高信息传播和存储的效率。在这里,我们考虑信号 x 是稀疏的,即 n 个维度中大部分元素为零,只有少量的非零元。
上面这个方程组 的目的就是利用压缩的信号 b,恢复原始信号 x。
如果你认为这就是一个简单的求解线性方程组问题的话,那就大错特错了,因为如前所述,这个方程组有无穷多个解!
实际上,原始信号重建问题对应的是一个约束问题:
即 在满足约束 的条件下,经可能地减少 x 中非零元的数目。
不幸的是,上述问题并不是一个凸优化问题,因为 表示非零元个数,不是一个凸函数。
取而代之,我们将优化问题中的目标函数换成
范数和
范数来看看,即考虑优化问题:
上述问题可以用现有的凸优化求解器快速求解! 因为
范数是凸函数,而替代问题2 可以通过一些变换转换成凸问题。
结果如下:
- 图1 对应原始的稀疏信号;
- 图2 对应在 范数约束下重建的信号;
- 图2 对应在 范数约束下重建的信号。
从结果可以看出, 正则化不能保证稀疏性,而 正则化可以!
以下是在 matlab 中调用 CVX 的 mosek 求解器,对上述 约束问题求解的代码。
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));