压缩感知与稀疏模型——Sparse Signal Model

这篇博客介绍第二节课的一些内容。虽然第二章题目是Sparse Signal Model,但是这篇博客还介绍了很多高维下的内容,因此内容是比较杂的。不过这个也就是上课内容的记录,而不是完全按照博客题目分的。

从马毅老师的课上目前学到的最重要的收获,就是遇到一个要优化的问题时从几何的角度来思考,这样能够辅助理解很多东西。

L1-norm为什么能得到更稀疏的点相对于L2来说?

之前的利用sparse representation来进行超分辨的文章中,有一个点是,我们希望得到最稀疏的解,但是L0范数的问题是NP-hard的,因此无法解决,所以将L0问题放宽到了L1问题。我们都知道L1范数能得到比L2范数更稀疏的解,为什么呢?

假如有下面待解决的问题:

对于0范数优化目标是最稀疏的解,是很容易理解的。因为0范数就是不为零的个数。通常解决方法是将上面的问题拓展到1范数:

现在我们想想,1范数的定义是什么:$Vert x Vert_1 = \Vert x_1 +…+x_n \Vert$。现在考虑三维的情况(更高的维度是类似的,但是我们画不出来),三维坐标下,L1范数的等高线是正8面体,而L2范数的等高线是球体的形状。而$y=Ax$,映射到三维坐标下则是一个平面。

我们在找满足$y=Ax$的情况下最小的$\Vert x \Vert_1$值,实际上就是这个等高线不断扩充,最早正八面体与平面接触的那个点。

而从概率的角度上,正八面体的6个顶点与平面接触的概率是最高的,而这几个顶点相对于其他的点来说是更为稀疏的(只有一个entry非零,如果是边的话则是有两个entry非零)。而球体(L2-norm)边界与平面接触概率是一样的,因此L1范数应该能比L2得到更稀疏的解。

$l^0,l^1,…l^\infty$的可视化

高纬度下的奇怪现象

我们一直想当然的任务高纬度下的数据与低纬度下的是一样的规律,但是实际上维度变高的情况下会出现一些反直觉的现象。

  1. 想想二维的圆和三维的球,对于更高的维度,也有一个类似这样的球。我们称围绕直径一圈的圆为赤道,赤道周围计算得到的体积,就占据了整个球的大多数(99%)。这个现象非常的反直觉,而且任意一个这样的赤道,都能得到相同的结论。

  2. 在维度变高情况下,$l^1$小于某个值相对于$l^\infty$占得体积越来越小,当维度比较高的时候接近于0。这个是比较好理解的。

还有一些想不起来了。总之在高维情况下,反直觉的现象会出现很多,因此很多的定理都是实践发现后才开始证明。这体现了编程验证能力的重要性。

L0 norm的sparse solution

下面考虑这样一个问题。加入我们有观测量$y \in \mathbb{R}^m$,并且知道矩阵$A$,有$y = Ax_0$,我们的目标是恢复$x_0$。如果我们知道$x_0$是稀疏的,那么选择最稀疏的解$y = Ax$是合理的,也就是我们需要解决的问题为:

这里定义一个符号:

一个向量的support是一个集合,包含了entry不为0的索引。那么上面的问题就是让我们找到一个有最小support的向量。最容易想到的方法就是尝试所有的可能性。support集合一定是集合$\{1,…,n\}$的子集,因此我们遍历集合$\{1,…,n\}$的所有子集,现在有support集合$I \subseteq \{1,…,n\}$,那么我们可以得到下面一组等式:

这里$A_I \in \mathbb R^{m \times \vert I\vert}$是矩阵$A$的子矩阵,由entry索引为集合$I$中的向量组成。这样我们尝试解出$x_I \in \mathbb R^{\vert I\vert}$。如果解存在,将其他的entry设为0即可。算法描述如下图:

利用这个算法做了随机实验,首先随机生成$x_0,A$,其中$x_0$是稀疏的。然后根据$y=Ax_0$得到观察量$y$,接着利用上面遍历的算法对$x_0$进行恢复。得到下面的图:

可以发现只要原先向量中非零个数$k$不要太大,也就是足够稀疏,基本上总是能成功找到解。上述实验的size是$5 \times 12$。出现这种现象的原因是什么,有没有数学上的解释?

为了解释为什么L0范数最小化能成功得到解,首先思考一下什么情况下是无法得到解的。假如$x_0 \in null(A)$,也就是属于$A$矩阵的零空间,那么:

对于上面的情况,最小化L0范数得到的解永远是0向量,无法恢复得到$x_0$。因此我们可以推测,如果$null(A)$包含了稀疏向量,那么我们可能无法恢复原来的稀疏向量。实际上,反过来也是正确的,如果$null(A)$不包含稀疏向量,通过解决L0最小化总是能恢复足够稀疏的原向量。上面的陈述有点让人看不懂,用数学语言说明并证明一下:

现在有原向量$\Vert x_0\Vert \leq k$,并且有:

那么如果$x_0$就一定可以被恢复。

首先,假如我们估计量为$\hat x$,那么由于$\hat x$是根据最小化L0范数得到的,因此:

定义估计误差为:

可以得到:

而:

那么根据$(*)$,我们可以得到:$\epsilon = \sigma = 0$。也就是最小化L0范数成功恢复了原来的稀疏向量。

而$(*)$是矩阵$A$的一个性质,因此一个好的矩阵$A$应该是零空间没有稀疏向量的。实际上,要保持这个性质的充要条件是,任意矩阵A的$2k$列是线性独立的。

对于上述现象在线性代数中有一个定理,引入矩阵新的概念Kruskal rank($krank(A)$),指的是最大的$r$,$A$的任意$r$列都是线性独立的。

上述理论描述为:

假定$y = Ax_0$,而且有$\Vert x_0\Vert_0 \leq \frac{1}{2} krank (A)$,那么$x_0$也是唯一下面最小化问题的唯一最优解:

Projected Subgradient Decent for $l^1$ minimization

对于$l^1$的最小化的解决实际上与我们经常遇到的$l^2$范数是有较大区别的。现在有下面这样的最小化问题:

上面的问题有两个需要解决的困难:

  1. nontrivial constraints(非平凡约束):需要满足$Ax = y$.
  2. nondifferentiable objective(不可导的目标函数)。$l^1$问题不是可导函数,局部会出现不可导的情况。下面对这两个问题一一进行分析解决。

Constraints

解决Constraints的最简单的方法是使用投影梯度。所有满足约束的向量,构成一个子空间(subspace),我们计算出来的梯度可以投影到这个子空间上,利用投影之后的梯度来进行梯度下降,就是投影梯度下降的思想。一般来说,带有约束的优化问题可以写成下面的形式:

其中$C$是一个约束集合,一般的梯度下降算法迭代步骤为:

而投影梯度算法就是将结果$x_{k+1}$投影到$C$上。一个点$z$在集合$C$上的投影也就是集合上距离$z$最近的点,定义为:

对于一般的$C$,投影可能不存在或者不唯一,但是对于凸集合来说,投影是很好的被定义了,而且有很多有用的性质。如果$A$是行满秩的,那么集合$C = \{x|Ax = y\}$上的投影有一个非常简单的形式:

有兴趣了可以尝试推导一下上式的计算过程。

Nondifferentiability

对于不可导问题,因为我们优化的问题虽然不可导,但是依然是凸函数,因此我们可以使用次梯度。这里简单说一下什么是次梯度,实际上一张图就可以看明白了:

从几何的角度上来说,我们希望梯度提供一个迭代的方向,但是由于不可导,梯度并不存在,但是实际上很多方向都可以作为迭代方向的估计,如上图中,红线的斜率就是次梯度。

对于凸函数,次梯度$c$需要满足:

可以看到次梯度一般不是唯一的而是一个区间,实际上我们对不可导的点左右求极限:

那么处于$[a,b]$之间的都可以称为次梯度。

对于次梯度更严格的定义如下:
$f: \mathbb R ^n \rightarrow \mathbb R$是一个凸函数。那么$f$在$x_0$点的次梯度为满足下式的任意$u$:

而次导数是在该点所有次梯度的集合:

因此对于次梯度下降,迭代步骤为:

截止到这里,我们就可以推导出次梯度下降算法了:

由于次梯度是无穷多个的,对于一维的变量,$\vert x\vert$在$x=0$的次梯度为区间$[-1,1]$的任何一个数字都可以,也就是:

我们用$\partial \vert \cdot \vert(x) $来表示次导数。对于多维向量,有下面的引理:
Lemma: 现在有$x\in \mathbb R^n$, $I = \text{supp}(x)$,那么:

这里,$\boldsymbol{P}_{I} \in \mathbb{R}^{n \times n}$是坐标$I$的正交投影:

证明这里就不赘述了。

到了这里我们可以总结出次梯度下降算法,用来解决$l^1$的minimization(注:下图中$A^* = A^T$,表示转置):

Exprimental Result

都说用$l^1$来解决稀疏解是相对有效的方法,因此在这里类似于$l^0$,做了一个实验,看$l^1$对恢复稀疏向量上表现如何。我们知道在$l^0$的情况下,只要满足一定条件它一定能恢复原来的稀疏向量。使用同样的方法设计实验,维度是$200\times 500$,非零entry个数$k$从0到500,算法使用上述提到的投影次梯度下降。这其中,首先随机生成$A,x$,根据$y=Ax$得到$y$,再试图用上述算法去恢复$x$。由于数值精度问题,计算机中完全等于0是比较困难一件事,因此设定一个阈值,只要小于该阈值就视为零。算法的代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
# A implementation of Projected Subgradient Decent
# created by Guoqing Zhang 2019.7.14

import numpy as np
from numpy.linalg import inv
import datetime
import matplotlib.pyplot as plt
# Use T and tildeX to compute projection
def computeProjection(A,y):
I = np.identity(A.shape[1])
I = np.mat(I)
A = np.mat(A)
y = np.mat(y)
temp = A.transpose() *inv(A * A.transpose())
T = I -temp*A
tildeX = temp*y
return [T,tildeX]
def randomSequenceNoRepeat(k,n):
r = list(range(0,n))

for i in range(k):
j = np.random.random_integers(0,n-1)
#print(j)
r[i],r[j] = r[j],r[i]
return r[:k]
# k means the sparsity, m, n is the dimension of matrix A
def randomGeneratexAy(k,m,n):
A = np.random.randn(m,n)
x = np.random.rand(n)
A = np.mat(A)
x = np.mat(x).transpose()
zero_indice = randomSequenceNoRepeat(n-k,n)
for i in zero_indice:
x[i] = 0
y = A*x
return x,A,y
# recoverX
def recoverX(A,y):
T,tildeX = computeProjection(A,y)
x = np.zeros(A.shape[1])
x = np.mat(x).transpose()

for i in range(1,10000):
#print(i,"iteration...")
x = tildeX + T*(x - 1/i * np.sign(x))
return x
def drawResult(sf):
x = list(range(0,len(sf)))
plt.plot(x,sf)
plt.show()

if __name__ == "__main__":
m = 200
n = 500
frac = [0 for i in range(501)]
last = 500
for i in range(50):
frac[i] = 1.0

for k in range(50,501):
success = 0
for i in range(60):
#starttime = datetime.datetime.now()
x,A,y = randomGeneratexAy(k,m,n)
#print(x)
hatX = recoverX(A,y)
hatX = hatX - x
hatX = np.fabs(hatX)
hatX[hatX<0.001] = 0
hatX[hatX!=0] = 1
#endtime = datetime.datetime.now()
#print((endtime - starttime).microseconds/1000.0)
if np.sum(hatX) == 0.0:
success = success + 1
frac[k] = success/60
print(k,success/60)
if(frac[k] == 0.0):
last = k
break
for i in range(last,501):
frac[k] = 0.0
frac[0] = 0.0
drawResult(frac)

最后结果如图:

这个图不是非常的顺滑,不过有了大致的趋势,可能是迭代次数不够导致的,也不知道是不是实现的问题。一般来说这个曲线如下图:

可以看到在稀疏性很大时候,$l^1$也几乎总是能恢复到原始的向量。这体现了$l^1$范数在解决稀疏问题的优越性。不过上述程序运行速度很慢,因为投影次梯度收敛速度很慢。之后会介绍更快速的优化算法。