分布式压缩感知的基础向内容梳理

写在前面的一点废话

    
2020年4月27日的一点感想：
    在经历了传统数据稀疏这个硬骨头之后，我发现我自己太想当然了，有些事情想的太过于简单了，我要做的事情是我自己选的，真的要摆正心态，不能着急！本人才疏学浅，能力有限，有误之处还望不吝指教。
    
2020年5月1日的一点感想：
    感觉自己就像是一个有坏道的硬盘，跑不起来最高的速度，有些东西存进去了就拿不出来，一点一点的修复，一点一点的补上坏的地方。
    

OMP与SOMP内容梳理

目前遇到的一些问题

    目前我做了对原始信号的稀疏度分析、对小波变换后的信号做了观测和重构，信号本身的特性发掘不出来，多路重构仅在重构所需数据量上与单路有所下降，并无大的创新。完成了OMP单路算法的重构；DCS-SOMP多路联合重构的完整过程。整个实验中缺少其他算法的对照组；缺少创新性；对于符合原始信号特性的联合稀疏不明确。
    以下，准备从寻找信号特性的联合稀疏和联合重构算法两个方面入手，首先先梳理一下现在的一些算法流程，巩固一下，也便于后续文字叙述直接使用。

明确几个概念

何为内积？

    此处的内积是指向量之间的数量积，几何意义是一个向量投影在另一个向量上的大小，在代数计算上就是向量A=（x₁，y₁），向量B=（x₂，y₂），A与B的内积写做A·B=x₁*x₂+y₁*y₂，若引申到N维则是x₁x₂…x_n+y₁*y₂…y_n,多个向量相互做内积也是如此。
    几何空间上A·B=|A||B|cosθ.
    内积可以表示两个向量之间的相关趋势，两个向量垂直则内积为0，两个向量平行内积为1,所以在迭代过程中通过传感矩阵与各列残差的内积来选出一个最大值，作为最相关的依据，此处后面详细说明。

何为投影？

     从二维角度上来说，我们常见的投影就是
当然一般的投影是位于坐标轴上的，在此处我们要说明的是一个向量 $b$ 在向量 $a$ 上的投影向量 $p$。

显然 $p$ 在 $a$ 上，所以 $p=xa$ ，$p$ 为 $x$ 倍的 $a$ $b$点与$p$点之间会存在一个误差$e=b-p$，且向量$e$是垂直于$a$向量的，由$p$是$a$的列空间中的一个子空间，所以垂直于列空间的是左零空间，所以$a^{T}e=0$成立，展开就是 $$a^T(b-xa)=0\Rightarrow a^{T}b-a^{T}xa=0\Rightarrow x=\frac{a^{T}b}{a^{T}a}$$
代入到 $p=xa$ 中，$$p=a\frac{a^{T}b}{a^{T}a}\Rightarrow p=\frac{a{a}^T}{a^{T}a}b$$
则可得到$p=Pb$，也就是说$b$乘了一个矩阵就变成了他的投影，所以我们称$P$为投影矩阵。
$$P=\frac{a{a}^T}{a^{T}a}$$

     当整个维度上升到2维以上时，就是说一个向量$b$，投影于A矩阵各列张成的列空间中，对于此处参见线性代数数学基础中对于四个基本空间部分的论述。此时的投影矩阵是$$P=A(A^{T}A{)}^{-1}{A}^T$$一般来说A不一定是方阵，所以其也不可逆，故$A^{T}A$为方阵才存在可不可逆的说法，这个式子中间也不可展开。

何为最小二乘法?

    通俗但不严谨的讲，最小二乘法是一种拟合散点使生成的曲线与各个数据点之间误差的平法和最小的方法。是一种将非线性数据按线性拟合的方法 （此句可能有歧义）。
    举个例子说，现在在平面直角坐标系上有4个点，A(1,5)、B(2,7)、C(4,8)、D(5,10)，想找一种方法来表示这4个点的变化趋势，此时我们设定y_i=ax_i+b，将4个点分别代入
                            5=a+b；7=2a+b；8=4a+b；10=5a+b；
此时我们想求得一组a和b使得上述4个式子均成立，但显然是不可能的，所以我们换种思路，使得没个点到拟合曲线的误差最小，设误差为   ε=Σ(y(x_i)-y_i(x))²=(a+b-5)²+(a2+b-7)²+(a4+b-8)²+(a*5+b-10)²
此时获得了所有点与目标曲线之间的一个差值的平方，代表了每一个散点与目标曲线之间的绝对距离，如果这个距离最小，则说明我们的曲线到每一个点的距离都是最小的，即便不能穿过，但是可以依照整体的趋势划分。整理后得到了ε(a,b) 的表达式，对该式子，我们分别求取a，b的偏导数，导数表示一个函数的变化率，偏导数是当这个函数不止有一个自变量时每个自变量的变化趋势，所以如果想取误差的最小值，我们需要分别求取当ε关于a和关于b的偏导数，当一阶导数出现0时说明函数出现拐点，为极值点（此处仅考虑示例中的情况为极小值点）
可列出两个式子：
$$\frac{\partial }{\partial a}\epsilon (a,b)\frac{\partial }{\partial b}\epsilon (a,b)$$
此时可以作为正常方程求解a和b，求出a=a₁，b=b₁即可求的一个式子y=a₁x+b₁，这个式子与各点的误差平方和最小，所以依据这种标准求取曲线的方法叫做最小二乘法。
    以上叙述的是最小二乘法在微积分代数领域中解超定方程组的一个示例，个人水平着实有限。。。不求严谨，但求叙述通俗。
    实际上我们在压缩感知中是要解欠定性方程组，同样的，可以用最小二乘法去拟合一个曲线，实现与目标点值的逼近，此处也涉及投影矩阵的概念。
    
    首先，我们确定一个要求解的问题是
$$\{\begin{array}{rlr}{y}_1& =& a_{1,1}{\theta }_1+a_{1,2}{\theta }_2+a_{1,3}{\theta }_3\\ y_2& =& a_{2,1}{\theta }_1+a_{2,2}{\theta }_2+a_{2,3}{\theta }_3\end{array}$$
此方程由两个方程联立，具有3个未知数，系数矩阵为
$$A=\left(\begin{array}{c}{a}_{1,1}，a_{1,2}，a_{1,3}\\ a_{2,1}，a_{2,2}，a_{2,3}\end{array}\right)$$
即 $y$ = $A$ $\theta$，使用最小二乘法求解该方程就是
$$\hat{\theta }=(A^{T}A{)}^{-1}{A}^{T}y$$
最终可以获得一个近似的 $\hat{\theta }$ 作为最优解即为所求，这个是在线性代数上做最小二乘。两者的联系见线性代数基础部分详解。

说明一下本文的符号体系

    此处是以压缩感知的常规形式给出，$x$_N*1 为原始信号，$\theta$_N*1为稀疏信号，$y$_M*1为观测信号，$\varphi$_M*N为观测矩阵，$\psi$_N*N稀疏矩阵，K为信号稀疏度，N为原始信号长度，M为观测信号长度，A为传感矩阵。
                                    $y$_M*1=$\varphi$_M*N $\theta$_N*1   A=$\varphi$_M*N $\psi$_N*N

其中N>>M>K，$\psi$_N*N一般为正交矩阵，也就是$\psi$^H_N*N = $\psi$^-1_N*N
后文中出现一些新符号时可能会就地解释。。。完善的时候再修改此处

OMP算法流程

    OMP的意思是Orthogonal Matching Pursuit，正交匹配追踪。
    关于他的前世今生不多说，是一种经典算法，前身是MP算法，加入了正交的步骤是使得求解内积时为0，不再重复计算残差，加快了算法的收敛速度。
    OMP的算法输入部分：观测信号y，传感矩阵A，信号稀疏度K。
    OMP的算法输出部分：重构的估计稀疏向量 $\hat{\theta }$ ,残差值 $r$_N*1。（这个$r$中只含有k次迭代的残差）
    首先说明，OMP算法并不是专门应用于压缩感知重构信号的，他是一种用于解欠定性方程组的算法，其他的重构算法也都是在解决欠定性方程组的问题，如何在欠定性方程组中找到一个最优解，也就是l0范数的问题，此处我们先不谈l0范数的问题，单纯从算法角度出发，结合一些前述知识去看他的运行步骤和算法流程。

1.  输入变量$y$，$A$，K，已知 $y$_M*1=$A$_M*N $\theta$_N*1 ,我们要求解的问题就是在已知 $y$ 和 $A$ 的基础上求 $\theta$，显然这是一个求解欠定方程的问题。因为$\theta$的维度大于$y$的维度。设定一个残差 $r$ 作为求解的目标，迭代之初，我们与目标函数的残差就是 $y$ 本身，即 $r$₀=$y$ ( $r$₀的角标表示迭代次数 )要在迭代的过程中使残差最小，也就让我们拟合的 $\hat{\theta }$ 更加接近于真实的 $\theta$。还需初始化一个序号索引集
                 （---------------------称该步骤为算法初始化过程---------------------）

2.  取传感矩阵 $A$ 的每一列（共N列）分别与残差 (也是个列向量) 做内积，可以获得一个N维的向量，在此向量内寻找一个绝对值的大的项，也就是投影值最大的项，认为该列与目前的残差最为相关(属于个人理解，可能不严谨)，记录下其位置序号以及 $A$ 中该序号列的数值。

3.  将这一列数记做 $a$_Max 作为方程 $y$_M*1=（$a$_Max）_M*1 $\theta$_N*1 中的系数，此时该式为一个欠定性方程组（通俗理解就是解的个数大于所联立方程的个数）此处不好理解可将$a$_Max这一列向量与$\theta$这一列向量展开，可以得到$$\left(\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_M\end{array}\right)=\left(\begin{array}{c}{a}_{1,Max}\\ a_{2,Max}\\ .\\ .\\ .\\ a_{M,Max}\end{array}\right)\ast \left(\begin{array}{c}{\theta }_1\\ {\theta }_2\\ .\\ {\theta }_M\\ .\\ {\theta }_N\end{array}\right)=\{\begin{array}{rlr}{y}_1& =& a_{1,Max}{\theta }_1+a_{1,Max}{\theta }_2+...+a_{1,Max}{\theta }_N\\ y_2& =& a_{2,Max}{\theta }_1+a_{2,Max}{\theta }_2+...+a_{2,Max}{\theta }_N\\ .& & \\ .& & \\ y_M& =& a_{M,Max}{\theta }_1+a_{M,Max}{\theta }_2+...+a_{M,Max}{\theta }_N\end{array}$$
也就变成了$$\left(\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_M\end{array}\right)=\left(\begin{array}{c}{a}_{1,Max}\\ a_{2,Max}\\ .\\ .\\ .\\ a_{M,Max}\end{array}\right){\theta }_1+\left(\begin{array}{c}{a}_{1,Max}\\ a_{2,Max}\\ .\\ .\\ .\\ a_{M,Max}\end{array}\right){\theta }_2+...\left(\begin{array}{c}{a}_{1,Max}\\ a_{2,Max}\\ .\\ .\\ .\\ a_{M,Max}\end{array}\right){\theta }_N$$
通过前述所说的最小二乘法求解($a$_Max^T.$a$_Max)^-1$a$_Max^T$y$= $\hat{\theta }$。

4.  该 $\hat{\theta }$为这一次迭代过程的结果，我们认为这个$\hat{\theta }$是当前的一个最优解，并用他重新求取
$\hat{y}=A_{索引列的集合}\hat{\theta }$ $就是A_{\Omega }$，随着迭代次数增加，$A_{\Omega }$的列数在不断扩大。将$\hat{y}$ 与原值$y$做差，作为更新残差 （$r$₁）=$y-\hat{y}$。此处可以判断残差是否符合要求精度，如果不符合将该残差继续代入步骤2中，迭代次数是输入的稀疏度，最多重构稀疏度K个。
（PS：但此时随着迭代次数的增加，$a$_Max的列数在不断增多，每迭代一次，就增加一列，相当于在保留了原有趋势的基础上加入修正余项，$\theta$的行数也随着迭代次数的增加不断增加，最终$a$_Max会变成一个M行K列的矩阵，依旧是一个扁矩阵，前述说过M>K，$\hat{\theta }$会变成一个K行的列向量，这个地方是我本人最开始理解的不好的地方，主要是数学底子太薄弱。。。）
重复 2、3、4 过程，不断更新残差，最终获得K个残差值。
                 （---------------------称以上步骤为算法贪婪迭代过程过程---------------------）

5.  最终获得了K个列数索引，而原 $\theta$ 是一个N长的向量，只不过是K稀疏的。所以我们按照这K个索引将最终所求得的 $\hat{\theta }$_k 排好顺序填进所有的非0点中作为最终重构的稀疏向量 $\hat{\theta }$输出。
                 （---------------------称该步骤为算法输出最终结果---------------------）
 

SOMP算法流程

     SOMP的意思是Simultaneous Orthogonal Matching Pursuit 同步/同时正交匹配追踪算法 该算法属于OMP拓展到多信号的一种方法。当信号数为1时，SOMP与OMP的步骤相同，此方法出现于

输入部分：J个观测信号的矩阵$S_{M\ast J}$，公共稀疏度K，传感矩阵$A_{M\ast N}$。
输出部分：重构的估计稀疏向量 ${\hat{\theta }}_{J\ast N}$ ,残差值 $r$_N*J。

1.  输入变量$S$，$A$，K，初始化残差为 $R_0$ = $S$ ,初始化索引集 $\Omega =\varphi$,初始化迭代计数变量i=1。
                 （---------------------称该步骤为算法初始化过程---------------------）

2.  通过解决一个简单的优化问题来寻找支撑集
$$\underset{\omega \in N}{\mathrm{argmax}}\sum _{J=1}^J|<r_{i-1}{e}_{\omega },A_{\omega }>|$$
     求解关于所有信号的残差与传感矩阵之间内积的和，就是每一个信号分别求内积，将多个信号的内积绝对值加和取最大值，作为本次迭代最相关列索引的坐标。并将该系数加入至索引集 $\Omega$ 内。
     此处与OMP的区别就在于OMP是输入信号是1列，将其与传感矩阵的各列做内积，取最大值，SOMP中输入的是个矩阵，将S矩阵的每一列分别与传感矩阵的一列做内积，加和，再分别与传感矩阵的下一列做内积，再加和，依次往复。
     式中的$e_{\omega }$为第$\omega$列的正交基。可以直接忽略他的作用，是为了表示在当前索引中张成的向量空间。

3.  求取投影，原文中给的说明就是求取一个投影系数，实际上就是和OMP中最小二乘法求取投影的步骤相同，对于每一个信号在该列上求取一个近似解。
每一个信号都会获得一个$\theta$故为${\hat{\theta }}_{i\ast J}$（i为迭代次数）

4.  分别更新每一列的残差值 $r_{N\ast J}$=$S_{M\ast J}-A_{\Omega }$* ${\hat{\theta }}_{i\ast J}$ ，对残差的大小进行判断，是否满足重构精度要求，不满足可重复 2、3、4 过程，不断更新残差。若满足了可以终止迭代。最终至多获得K个残差值。
                 （---------------------称以上步骤为算法贪婪迭代过程过程---------------------）

5.  最终获得了K个列数索引，而原 $\theta$ 是一个N长的向量，只不过是K稀疏的。所以我们按照这K个索引将最终所求得的 $\hat{\theta }$_k 排好顺序填进所有的非0点中作为最终重构的稀疏向量 $\hat{\theta }$输出。
                 （---------------------称该步骤为算法输出最终结果---------------------）

总结SOMP

     SOMP的重构问题在于重构多个信号的时候要确保所有信号的稀疏位置相同，因为会求取多次最小二乘系数，所以各个信号的幅值不同没有关系，各次运算之间是独立的。如果出现信号与信号之间稀疏位置不同的情况，在计算支撑集时就会产生支撑集选择错误的情况，造成重构不准确。
例如：6个信号联合重构时在该位置上有一个4个信号有值，2个信号的值为0，则计算内积时只会计算4个信号与传感矩阵的内积，会增加序号索引不正确的可能性。
     但是实际仿真中发现这种情况的重构效果还可以，在2个信号的0位置上会出现一个较小的值，但是与整体幅度趋势区别较大，当信号整体稀疏系数（稀疏信号的幅值）较小且观测数（M）不足、稀疏位置不同的情况下容易出现索引位置不正确的情况。
     SOMP的一个好处是引入了多信号之间在公共稀疏基下的稀疏位置相同这一相关性，仅用了内积绝对值加和这一方法就在一定程度上加强相关度最大列的索引依据，使得联合重构时即使观测点数较少，也可以重构出原始信号，这一点在OMP上有了很大改善，OMP中的M经验值要为3-4倍的K效果较好，而SOMP中对于联合多路信号，信号数越多，我们需要的观测数就可以越少，小于2都可以，本人在实际试验中采用M=1.5倍的K进行6路信号重构时，误差就可以控制在2%以内，相比于OMP好了很多，在DCS的论文中，给出了一个趋势就是当信号数趋于无穷时，且各个信号在公共稀疏集上稀疏位置相同，则观测数可趋近与M=K+1即可。

DCS-SOMP算法流程

From【BARON D, WAKIN M B, DUARTE MF, et al. Dis-tributed compressed sensing[J]. Neural Information Processing Systems(NIPS),2005,22(10):2729-2732.】
    DCS-SOMP与SOMP有所不同，文中说DCS-SOMP是对SOMP在多信号拓展时做了适应性调整，这个说法我没有在文中找到实际的依据，但是不同的是DCS-SOMP对正交投影的运用方法与OMP和SOMP有些写法上的区别，其实方法本身应该是相似的，只是表现方法不同，在DCS-SOMP中加入了专门的正交化步骤和去正交化步骤。
这个步骤目前我梳理清楚了，但是其中涉及施密特正交化与最小二乘法之间的联系，正在补充线性代数的相关知识，后续更新这里，争取讲明白。
    
                    叙述过程中部分对原文符号体系进行了修改，方便与前文对齐。

输入部分：J个观测信号的矩阵$S_{M\ast J}$，公共稀疏度K，传感矩阵$A_{M\ast N}$。
（NOTES： 原文中的输入时观测矩阵，其默认进来的向量就是稀疏的，我们此处可以想象原始信号就是一个变换域的稀疏信号，且文中说明了每一列信号J的随机观测矩阵都是不同的，由于每一个观测矩阵都是随机矩阵。且稀疏矩阵与随机矩阵组成的传感矩阵皆默认满足RIP准则，此处为了说明方便讲所有信号采用同一个随机观测矩阵进行观测。）
输出部分：重构的估计稀疏向量 ${\hat{\theta }}_{J\ast N}$ ,残差值 $r$_N*J。

1.  输入变量$S$，$A$，$K$，初始化残差为 $R_0$ = $S$ ,初始化索引集 $\Omega =\varphi$,初始化迭代计数变量i=1。
                 （---------------------称该步骤为算法初始化过程---------------------）

2.  通过一个优化问题 $$\underset{\omega \in N}{\mathrm{argmax}}\sum _{j=1}^J\frac{|<r_{j,i-1},A_{\omega }>|}{||A_{\omega }|{|}_2}$$                                   式中的$\omega$表示观测矩阵的列数

     这个优化问题其实是求取每个信号的残差在观测矩阵中每一列的投影，参考前文说的投影部分，其实本意也是去最相关的那一列向量作为后面求解的支撑集，将该列的序号加入 $\Omega$。

3.  将索引到的那一列数据进行正交变换
$${\gamma }_i=A_{j,n_i}-\sum _{t=0}^{i-1}\frac{⟨A_{j,n_i},{\gamma }_t⟩}{{\Vert {\gamma }_t\Vert }_2^2}{\gamma }_t$$随着迭代次数的增多，新加入支撑集的索引列数据依次正交化，并将正交后的数据保存为 ${\gamma }_i$，此处$i$为其中$i$为$\gamma$的第$i$列，$i$表示迭代次数，角标$t$也同理，表示在施密特正交化过程中 $\gamma$ 索引列的序号，文中默认每一个信号的观测矩阵不同，所以给出了每一个信号不同的正交矩阵，此处我们采用的观测矩阵为相同的，故只产生一个正交矩阵，或者说每一个信号到这步的正交矩阵都是相同的。

4.  开始更新求取的估计值 $\hat{\beta }$ 和残差 $r$ ，注意此处的$r$和$\gamma$不是一个东西。。。最开始看的时候一时粗心看成了一个元素，推算了好久也没搞明白什么意思，最后发现他俩是两个完全不同的符号，输出出来后$r$被拉长了，就长得较为相似，吃了粗心的亏。
$$\begin{array}{rl}{\hat{\beta }}_j(i)& =\frac{⟨r_{j,i-1},{\gamma }_i⟩}{{\Vert {\gamma }_i\Vert }_2^2}\end{array}$$

$$\begin{array}{rl}{r}_{j,i}& =r_{j,i-1}-\frac{⟨r_{j,i-1},{\gamma }_i⟩}{{\Vert {\gamma }_i\Vert }_2^2}{\gamma }_i\end{array}$$
     此处估计值 $\hat{\beta }$ 的更新策略是取残差在索引集正交变换后的矩阵上的投影，每一个信号与公共的索引向量单独求取估计值，故一共有$j$个估计值，且随着迭代次数的增多，${\hat{\beta }}_j$ 的维度也在扩大， 文中也称该估计值为正交系数向量。
     此处估计值 $r$ 的更新策略其实和OMP类似，OMP是取估计值 $\hat{\theta }$与 $A_{\Omega }$ 中现有的所有索引集做乘法，获得当前迭代过程中的近似解，DCS-SOMP中是取了 $\hat{\beta }$与 $A_{\Omega }$ 的索引集做正交变换后的 $\gamma$ 做乘法，与原残差做差，更新残差值。其实这一步和上文中说的OMP及SOMP的区别就是DCS-SOMP全部是在正交集上操作，本质原理上是一回事。
5.  检查残差阈值是否符合精度要求，如果符合则停止迭代，不符合则重复上述2、3、4过程，该算法至多迭代M次，即观测点数为最大迭代次数，（文中说这个M次是局限于第三步中的正交策略，正交变换只在每一次迭代中才会更新一列新的正交索引值，但是这类算法的最大迭代次数不应该都是稀疏度K么？即便是迭代了观测数M次，也应该只有K个有效值。此处存疑）最终获得一个 $\hat{\beta }$和 ${\gamma }_i$ 两个下一步去正交化的关键参数。
                 （---------------------称以上步骤为算法贪婪迭代过程过程---------------------）

6.  使用QR分解对$A_{\Omega }$进行分解 $A_{\Omega }=Q\ast R$ 矩阵的QR分解是将一个矩阵分成Q（正交矩阵）和R（上三角矩阵）两个矩阵相乘的形式，该Q矩阵前面的步骤中已经获得，就是 $\gamma$ 的集合，因此我们可以求得R矩阵。
     输入的观测信号 $y$=$\gamma$ 的集合*$\hat{\beta }$的集合 ，也就是 $y$=Q*$\hat{\beta }$，且 $y$=$A_{\Omega }$*${\hat{\theta }}_{\Omega }$，
代入后可得 $Q\ast \hat{\beta }$=$Q\ast R\ast \hat{\theta }$, Q消去，可求得
$${\hat{\theta }}_{\Omega }=R^{-1}\ast \hat{\beta }$$
最后，依旧按照索引序号将所求的 ${\hat{\theta }}_{\Omega }$分配至索引位置，输出最终的 $\hat{\theta }$。

                 （---------------------称该步骤为算法输出最终结果---------------------）

总结DCS-SOMP

     DCS中在索引到的序列中加入了正交变换的步骤，使得整个运算进行在A的列空间中，最终获得的也是投影系数，最后利用正交矩阵进行QR分解获得最终输出值。从本质上来看求解过程是一致的，就是多了正交化和去正交化的过程，整个过程用矩阵表示的更加宏观了。