正交匹配追踪

数字图像处理 专栏收录该内容
10 篇文章 0 订阅

这篇博文是在对Koredianto Usman《Introduction to Orthogonal Matching Pursuit》文章的翻译,后面附带了一些总结.

这篇文章是前面《Matching Pursuit (MP)》文章的继续. (注:原文中还是有一些小细节错误,请大家睁大眼睛阅读)

简介

考虑下面的问题:给定 x=[1.210] A=[0.7070.7070.80.601] ,计算 y=Ax .

这相当简单!

y=Ax=[0.7070.7070.80.601][1.210]=[1.650.25]

现在是比较困难的部分!给定 y=[1.650.25] A=[0.7070.7070.80.601] ,如何找到原始的 x (或者尽可能接近原始的x)?

在压缩感知(Compressive Sensing)术语中,从 x A中得到 y 叫做压缩,反之,从y A 得到 x 叫重建(个人感觉原文中好像说反了,还是我理解错了?). 重建问题不是一个很简单的问题。

一般地,x被称为original signal或者original vector, A 被称为compression matrix或者sensing matrix, y 称为compressed signal或者compressed vector.

基本概念

和前面在MP算法中讨论过的一样,感知矩阵A可以看做列向量的集合:

A=[0.7070.7070.80.601]=[b1b2b3]

其中,

b1=[0.7070.707],b2=[0.80.6],b1=[01]

这些列被称为基(basis)或者原子(atom).

现在,我们令 x=x1x2x3 ,则

Ax=[0.7070.707]x1+[0.80.6]x2+[01]x3=y
.

从该方程可以看出 y 是使用x中的稀疏对原子的线性组合. 实际上,我们知道 x1=1.2 x2=1 x3=0 . 换言之,对于得到 y 值,b1的贡献值为-1.2, b2 为1,而 b2 为0.

OMP算法和MP算法类似,都是从字典中找出哪一个原子对 y 值的贡献最大,接下来是哪个原子的贡献值大,以此类推. 我们现在知道这个过程需要N次迭代, N 是字典中原子的个数. 在该例子中,N等于3.

计算原子的贡献值

字典中有三个原子,

b1=[0.7070.707],b2=[0.80.6],b1=[01]

y=[1.650.25]

因此,各原子对 y 的贡献值为

<b1,y>=[0.7070.707]T[1.650.25]=0.7071.65+0.7070.25=1.34

<b2,y>=1.17

<b3,y>=0.25 <script type="math/tex" id="MathJax-Element-41"> = 0.25</script>

显然,原子 b1 y 值的贡献最大为-1.34(按绝对值计算).

在第一次迭代过程中,我们选择b1作为基,基的系数为-1.34.

当然,我们也可以使用矩阵的形式计算点积:

w=ATy=0.7070.800.7070.61[1.650.25]=1.341.170.25

原子和 y 的几何意义如图1:

原子和感知向量示意图

计算残差

现在已经选定了第一个基b1=[0.7070.707],贡献值系数为 λ1=1.34 . 将该部分贡献值从 y 中减去,残差为

r=yλ1b1=[1.650.25](1.34)[0.7070.707]=[0.70.7]

残差的几何意义如图2:

残差的几何意义

重复迭代

经过第一次迭代,我们选择出了基 b1 . 将选择出的基加入新的矩阵 Anew ,即:

Anew=[b1]=[0.7070.707]

贡献系数矩阵可以写成 xrec=1.3400

第一个元素是第一个贡献值-1.34,对应着第一个第一个基 b1 .

残差值为 r=[0.70.7]

现在需要从剩下的基 b2 b3 中选择对残差 r 贡献最大的.

w=[b2b3]Tr=[0.980.7]

所以, b2 的贡献值更大一些是0.98.

我们将选择的基 b2 添加到 Anew 得到

Anew=[b1b2]=[0.7070.7070.80.6]

接下来的步骤和MP算法稍微有些不同. 这里,我们需要计算 Anew y 的贡献(而不是b2 r 的贡献)

我们使用最小二乘法(Least Square Problem)解决该该贡献值问题:

如何求得λ1 λ2 的值使得 λ1[0.7070.707]+λ2[0.80.6] 最接近 y=[1.650.25]

表示成数学语言为:

minAnewλy2

在这个例子中,即

min[b1b2]=[0.7070.7070.80.6][λ1λ2][1.650.25]2

对于 minAnewλy2 的解可以直接使用公式:

λ=A+newy

其中, A+new Anew 的谓逆, A+new=(ATnewAnew)1ATnew

(注:可参看我的博文《Moore-Penrose广义逆矩阵》和《矛盾方程的最小二乘解法》)

在我们的例子中 A+new=[0.7070.7070.80.6]+=[0.60620.71430.80820.7143]

可以使用MATLAB中的内置函数pinv()进行伪逆的计算.

得到 A+new 后,可以计算 λ

λ=A+newy=[0.60620.71430.80820.7143][1.650.25]=[1.21]

得到 λ 以后,更新系数矩阵 xrec

xrec=1.210

同样的, λ xrec 第一个元素和第二个元素对应着选择出来的第一个和第二个基.

得到 Anew λ 以后,我们可以计算残差

r=yAnewλ=[1.650.25][0.7070.7070.80.6][1.21]=[00]

此时,残差值为 [00] ,我们可以停止计算或者继续下一步计算(这里停止,可以节省一些计算工作).

最后一次迭代

这一步不是必须的,因为残差已经完全消除了(很多实现OMP的软件都需要输入稀疏度 K 参数,这样经过K次迭代以后,无论残差大小都会停止迭代).

重新组织信号系数 xrec=1.210 ,这刚好是原始的系数.

其他例子

给定 x=0312 A=0.80.20.20.30.4110.30.10.40.40.8

求得 y=Ax=2.70.14.5

现在给定 A y ,使用OMP算法求解x.

  1. 我们有4个基(原子): b1=0.80.20.2b2=0.30.41b3=10.30.1b4=0.40.40.8

  2. 由于基向量的长度不是1,所以我们首先进行标准化.

    b1^=b1/b1=0.80.20.2/(0.8)2+(0.4)2+(0.2)2=0.94280.23570.2357

    b2^=b2/b2=0.26800.35780.8940

    b3^=b3/b3=0.95350.28600.0953

    b2^=b2/b2=0.40820.40820.8165

  3. 标准化的基向量对 y 的贡献w

    w=A^Ty=0.94280.23570.23570.26800.35780.98400.95350.28600.09530.40820.40820.8165T2.70.14.5=1.50854.78522.11674.7357

  4. 第二个基向量 b2 贡献值最大,所以将 b2 加入到 Anew

    Anew=b2=0.30.41

  5. 利用最小二乘法计算 xrec

    Lp=A+newy=[4.28]

    因为 Lp 对应着第二个基向量,所以 xrec=04.2800

  6. 接下来计算残差

    r=yAnewLp=2.70.14.50.30.414.28=1.4161.6120.22

  7. 接下来重复第3步,计算 b1^ b1^ b1^ r 的贡献.

    w=A^Tr=0.94280.23570.23570.95350.28600.09530.40820.40820.8165T1.4161.6120.22=0.90321.79021.4158

    选择第二个贡献最大的基 b3 ,其贡献值为1.7902.

  8. 将选择的 b3 加入到 Anew

    Anew=[b1b2]=0.30.4110.30.1

  9. 用最小二乘法计算 Lp=A+newy=[4.17021.7149]

    这个 Lp 对应着 b2 b3 ,因此 xrec=04.1721.71490

  10. 接着计算残差

    r=yAnewLp=2.70.14.50.30.4110.30.1[4.1721.7149]=0.2661.05360.5012

  11. 重复步骤2,进行最后一次迭代

  12. 计算 b1 b4 的贡献值 w=[b1^b2^]r=0.94280.23570.23570.40820.40820.8165T0.2661.05360.5012=[0.61720.7308]

    这里 b4 提供了最大贡献值0.7308.

  13. b4 加入 Anew=0.30.4110.30.10.40.40.8

  14. 利用最小二乘计算 Lp=A+newy=321 ,对应的 xrec=0312

  15. 接着计算残差

    r=yAnewLp=2.70.14.50.30.4110.30.10.40.40.8312=000

  16. 我们的迭代到此为止,因为此时残差已经为0,重建的 x xrec=0312,和原来的信号相同.

需要注意的问题

通过上面的迭代计算过程,我们应该注意如下几点:

  1. OMP中最大贡献值的计算需要对基向量进行标准化处理,不是由原始基得到的.
  2. 如果给定的基向量已经是单位向量,则不需要进行标准化.
  3. 基向量对 y 值的贡献的计算是通过标准化以后的基向量和残差的点积进行计算的.
  4. 在MP算法中,重建系数xrec的计算是通过基向量和残差的点积进行计算的,在OMP算法中, xrec 的计算是通过最小二乘法得到 Anew 相对于 y 的系数得到的,即Lp的计算. Lp 向量中的值直接对应于 xrec 的相应位置. Anew 通过每次对基向量的选择得到. 这个过程是需要时间的,因此OMP比MP慢. (不是应该OMP收敛的快吗?)
  5. 残差 r 的计算通过原始的y值和 Anew Lp 得到.
  6. 迭代的次数最多等于 A 矩阵的行数M,或者如果给定了稀疏度 K ,则迭代K次. 如果 K<M ,则已知的 K 可以加快计算结束,如果K未知,则迭代 M 次.

说明总结

在正交匹配追踪OMP中,残差总是与已经选择过的原子正交的。这意味着一个原子不会被选择两次,结果会在有限的几步收敛。

OMP算法 步骤描述:

输入:字典矩阵A,采样向量 y ,稀疏度k.

输出: x K稀疏逼近 x^ .

初始化:残差 r0=y ,索引集 Λ0= t=1 .

循环执行步骤1-5:

  1. 找出残差 r 和字典矩阵的列Ai积中最大值所对应的脚标 λ ,即 λt=argmaxi=1,,N|<rt1,Ai>| .
  2. 更新索引集 Λt=Λt1{λt} ,记录找到的字典矩阵中的重建原子集合 At=[At1,Aλt] .
  3. 由最小二乘得到 x^t=argmin|yAtx^| .
  4. 更新残差 rt=yAtx^t t=t+1 .
  5. 判断是否满足 t>K ,若满足,则迭代停止;若不满足,则继续执行步骤1.
相关推荐
©️2020 CSDN 皮肤主题: 成长之路 设计师:Amelia_0503 返回首页
实付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值