稀疏表示求解:OMP(The Orthogonal Matching Pursuit Algorithm)
文章目录
一、一些定义与原理
1. L 0 L_0 L0范数
定义为向量中非零元素的个数,用于描述向量的稀疏程度。
2. L 0 L_0 L0范数约束的稀疏表示问题描述(OMP求解目标)
问题描述如下:
min
x
∣
∣
x
∣
∣
0
,
s
.
t
.
A
x
=
b
\min_x{||x||_0}, s.t. Ax=b
xmin∣∣x∣∣0,s.t.Ax=b
对于稀疏表示这一特定应用,该公式可解释为找到最稀疏的
x
x
x(通过
L
0
L_0
L0范数作为罚项约束),使得
b
b
b能够在特定的矩阵
A
A
A下(稀疏表示任务中称
A
A
A为字典)被
x
x
x表示。我们的目的是求解稀疏系数
x
x
x。
3. L 0 L_0 L0范数约束的稀疏表示问题的求解是NP难的
引入一个概念:
矩阵
A
A
A的spark,即矩阵中最小的线性相关组中列向量的个数。
具体来说,如下矩阵的rank为4,而spark为3(即标红的三列是线性相关的,也是最小的线性相关列向量组)。
如果想要了解spark在压缩感知中的应用,可以去搜一下。
根据Michael Elad等人的证明:
当且仅当
s
p
a
r
k
(
A
)
>
2
k
spark(A)>2k
spark(A)>2k时,可以通过上述最小0范数优化问题得到稀疏系数
x
x
x的精确近似。我们可以直观想象一下spark的求解过程:
首先明确,
s
p
a
r
k
>
=
2
spark>=2
spark>=2,存在两个向量相等的情况时取等号。
进而,从
n
=
2
n=2
n=2开始搜索所有的列向量组合,直到第
n
=
N
n=N
n=N次时某一组合中的向量线性相关,则
s
p
a
r
k
=
N
spark=N
spark=N。
显然这种搜索所有可能性的方法由于矩阵
A
A
A庞大的列向量个数(过完备性),使得很难在多项式时间内得到问题的解。因此研究者们致力于如何优化最小0范数优化问题的求解。
具体来说,如果矩阵
A
A
A有2000列,
N
=
15
N=15
N=15的话,则需要大约
7.5
×
1
0
20
7.5 \times{10^{20}}
7.5×1020年来求解。
二、OMP算法
1. 引入残差(residual)
根据上述最小0范数优化问题描述,可以对
A
x
=
b
Ax=b
Ax=b限制不那么严格,引入如下误差:
r
k
=
b
−
A
x
k
r_k=b-Ax_k
rk=b−Axk
OMP的原理就是选择下一个非零值
x
x
x,以尽可能减少残差
r
k
r_k
rk使得
A
x
Ax
Ax更接近于
b
b
b。这就将最小0范数优化问题转化为了一个迭代求解问题。
2. 具体步骤
借用一下Elad大佬的👇
-
第一步(1&2),从矩阵 A A A中选择最佳列 a i 0 a_{i_0} ai0,使得这一列与一个标量 z z z相乘得到的结果与残差 r k − 1 r_{k-1} rk−1的 L 2 L_2 L2范数小于其他列。
求解细节如下:要想求得满足公式1的标量 z z z,首先对公式求导(极小值)得到: z o p t = a i T r k − 1 a i T a i = a i T r k − 1 z_{opt}=\frac{a_i^Tr_{k-1}}{a_i^Ta_i}=a_i^Tr_{k-1} zopt=aiTaiaiTrk−1=aiTrk−1然后再将 z o p t z_{opt} zopt代入公式1得到: E ( i ) = ∣ ∣ a i T r k − 1 ⋅ a i − r k − 1 ∣ ∣ 2 2 = ∣ ∣ r k − 1 ∣ ∣ 2 2 − ( a i T r k − 1 ) 2 E(i)=||a_i^Tr_{k-1}\cdot a_i-r_{k-1}||_2^2=||r_{k-1}||_2^2-(a_i^Tr_{k-1})^2 E(i)=∣∣aiTrk−1⋅ai−rk−1∣∣22=∣∣rk−1∣∣22−(aiTrk−1)2由此,最小化 E ( i ) E(i) E(i)的问题就转化为了最小化 ∣ a i T r k − 1 ∣ |a_i^Tr_{k-1}| ∣aiTrk−1∣。
最小化 ∣ a i T r k − 1 ∣ |a_i^Tr_{k-1}| ∣aiTrk−1∣就好求了呀👇
-
第二步(3),将第 i 0 i_0 i0列放入 S k S_k Sk中。
-
第三步(4&5),更新 x k x_k xk,与新的残差 r k r_k rk。
求解细节如下: S − k S-k S−k是从 A A A中之前每一轮最优的 a i a_i ai组成的,为下图标绿的部分,而 x x x中标绿的部分与之对应。
更新 x k x_k xk的步骤则可以简化成只更新 S k S_k Sk对应的部分,即绿色部分,一次简化计算。我们可以提取这些对应列与值,组成 A s A_s As。按照如下公式来更新 x k x_k xk:
x k = min x ∣ ∣ A s x − b ∣ ∣ 2 2 x_k=\min_x ||A_sx-b||_2^2 xk=xmin∣∣Asx−b∣∣22对上式求导得到 x k = ( A s T A s ) − 1 A s T b = A s + b x_k=(A_s^TA_s)^{-1}A_s^Tb=A_s^+b xk=(AsTAs)−1AsTb=As+b
-
重复上述步骤如果最终误差( r k r_k rk)小于设定阈值,则停止。
3. 算法命名原因
每一轮的更新是为了让更新的误差 r k r_k rk与矩阵 A A A正交(Orthogonal),即满足 A s T ( A s x k − b ) = − A s T r k = 0 A_s^T(A_sx_k-b)=-A_s^Tr_k=0 AsT(Asxk−b)=−AsTrk=0
三、OMP算法变体
OMP算法后续发展了一些变体,在运算复杂度和准确度之间存在trade-off
1. Least-Squares OMP (LS-OMP)
OMP算法通过判断残差判断某一轮最优的
a
i
a_i
ai,但这很有可能产生次优解。为了确保每一轮都可以选到最优解,LS-OMP舍弃残差,对逐个将每一列放入上一轮的
S
k
−
1
S_{k-1}
Sk−1中,并对比所有可能产生的矢量结果,最后进行选择,具体算法如下图所示。
LS-OMP的算法复杂度要高于OMP,但是效果更好。
2. Matching Pursuit (MP) 与 Weak MP
- MP算法:
为了简化OMP算法,MP算法对于第k次更新不是从零开始,而是利用第k-1次更新的结果。具体来说,由于已经得到了 x k − 1 = min x k − 1 ∣ ∣ A s k − 1 x k − 1 − b ∣ ∣ 2 2 x_{k-1}=\min_{x_{k-1}} ||A_{s_{k-1}}x_{k-1}-b||_2^2 xk−1=xk−1min∣∣Ask−1xk−1−b∣∣22在更新 x k x_k xk的时候可以将第一步修改为: x k = min z ∣ ∣ A s k − 1 x k − 1 + a i 0 z − b ∣ ∣ 2 2 x_k=\min_z||A_{s_{k-1}}x_{k-1}+a_{i_0}z-b||_2^2 xk=zmin∣∣Ask−1xk−1+ai0z−b∣∣22可以看到,上述公式等号右侧的三项中,第一项减第三项正好是上一轮的误差,所以上述公式可以转化为 x k ( i 0 ) = min z ∣ ∣ a i 0 z − r k − 1 ∣ ∣ 2 2 x_k(i_0)=\min_z||a_{i_0}z-r_{k-1}||_2^2 xk(i0)=zmin∣∣ai0z−rk−1∣∣22该式可以通过求导求得闭式解。因此,MP将OMP算法的第4步修改为: x k = x k − 1 x_k=x_{k-1} xk=xk−1 x k ( i 0 ) = x k ( i 0 ) + a i 0 T r k − 1 x_k(i_0)=x_k(i_0)+a_{i_0}^Tr_{k-1} xk(i0)=xk(i0)+ai0Trk−1 MP算法很有可能重复选择同一列 a i a_i ai所以它的准确度要低于OMP - Weak MP算法:
给定 ∣ a i T r k − 1 ∣ |a_i^Tr_{k-1}| ∣aiTrk−1∣的期望阈值( 0 − ∣ ∣ ∣ r k − 1 ∣ 2 0-|||r_{k-1}|_2 0−∣∣∣rk−1∣2)超过这个阈值t倍则停止(根据E(i)>=0)。