压缩感知和稀疏恢复基础

2023-08-08

字数统计: 10.7k | 阅读时长≈ 43 分钟

0. 文献研读中知识点补充

0.1 托普利兹(toeplitz)矩阵

(一) 基本定义

Toeplitz矩阵又叫做常对角矩阵(diagonal-constant matrix)，指矩阵中每条自左上至右下的斜线上之元素都为同一常数的矩阵。例如下面就是一个Toeplitz矩阵的例子：

$\left[\begin{array}{ccccc} a & b & c & d & e \\ f & a & b & c & d \\ g & f & a & b & c \\ h & g & f & a & b \\ i & h & g & f & a\\ \end{array}\right]$

其正式定义托普利兹矩阵_百度百科 (baidu.com) 为 —— 设$\boldsymbol T = [t_{ij}] \in \mathbb C^{n \times n}$，如果满足$t_{ij} = t_{j-i}$，即：

$\boldsymbol T = \left[\begin{array}{ccccc} t_0 & t_1 & t_2 & \cdots & t_{n-1} \\ t_{-1} & t_0 & t_1 & \cdots & t_{n-2} \\ t_{-2} & t_{-1} & t_0 & \cdots & t_{n-3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ t_{-n+1} & t_{-n+2} & t_{-n+3} & \cdots & t_{0} \end{array}\right]$

则称$\boldsymbol T$为托普利兹矩阵(Toeplitz matrix)。

(二) 对称Toeplitz和Hermitian Toeplitz矩阵

在信号处理领域中经常遇到一种特殊的T型矩阵，它除了具有一般T型矩阵的特点外，还是一个对称矩阵，即满足$t_{-i} = t_i$，形式如下：

$\boldsymbol T = \left[\begin{array}{ccccc} t_0 & t_1 & t_2 & \cdots & t_{n-1} \\ t_{1} & t_0 & t_1 & \cdots & t_{n-2} \\ t_{2} & t_{1} & t_0 & \cdots & t_{n-3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ t_{n-1} & t_{n-2} & t_{n-3} & \cdots & t_{0} \end{array}\right]$

在复数域上，若一个复Toeplitz矩阵的元素满足复共轭对称关系，即$t_{-i} = t_i^*$，形式如下：

$\boldsymbol T = \left[\begin{array}{ccccc} t_0 & t_1 & t_2 & \cdots & t_{n-1} \\ t_{1}^* & t_0 & t_1 & \cdots & t_{n-2} \\ t_{2}^* & t_{1}^* & t_0 & \cdots & t_{n-3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ t_{n-1}^* & t_{n-2}^* & t_{n-3}^* & \cdots & t_{0} \end{array}\right]$

则称之为Hermitian Toeplitz矩阵。

这种特殊的对称Toeplitz矩阵可仅由矩阵的第一行元素完全确定，因此常简记为$\boldsymbol T = \boldsymbol T(\boldsymbol t) = \mathrm{Toep}[t_0, t_1, \cdots, t_{n-1}]$。

0.2 封闭解、解析解、数值解释义

(一) 封闭解

首先，封闭解(closed-form expression)就是一个数学表达式；这种数学表达式包含有限个标准运算。

什么是有限？ 比如$x+y$就是一个有限表达式(只有两个变量、一个加法运算)；二次方程的根$\dfrac{-b \pm \sqrt{b^2-4ac}}{2a}$也是一个有限表达式；$\sum\limits_{i=1}^5 2^i$也是有限表达式(这里包含了5次求和运算)。但是$\sum\limits_{i=1}^{+\infty} 2^i$就不是一个有限表达式，因为他是无穷求和。

什么是标准运算？根据上面的定义，标准运算包含常量、变量、常见运算(如加、减、乘、除)以及一些函数(如$N$次根、指数、对数、三角函数、双曲函数)；但是极限、差分、积分都不能算作标准运算。

有了这两条标准，我们就能很容易地判断出一个解是不是闭合解了，如下示例：

$x+y$是一个闭合解，它仅包含变量和加法运算；
$\dfrac{-b \pm \sqrt{b^2-4ac}}{2a}$是一个闭合解，它包含变量、加、减、乘、除以及根号运算；
$\sum\limits_{i=1}^5 2^i$是一个闭合解，它包含变量、加法、除法以及指数运算；
$\sum\limits_{i=1}^{+\infty} 2^i$不是一个闭合解，因为它包含了无限次加法运算；
$\Phi(x) = \displaystyle\int_{-\infty}^x \dfrac{1}{2\pi} \exp\left(\dfrac{t^2}{2}\right) \text{ d}t$表示标准正态分布的累积分布函数，如果一个解中包含$\Phi(\cdot)$运算符，那么这个解就不是闭合解(积分运算不能算作标准运算)。

两点说明：

有些情况下，非闭合解可以转化为闭合解。比如无穷级数$\sum\limits_{i=1}^{+\infty} \dfrac{x}{2^i}$本身不是一个闭合解，但是利用等比数列的求和公式我们可以得到$\sum\limits_{i=1}^{+\infty} \dfrac{x}{2^i} = 2x$。因此，可以将非闭合解转化为闭合解。
维基百科上关于Closed-form expression的定义的最后一句提到：“The set of operations and functions may vary with author and context”。换句话说，关于标准运算的定义并不是非常的严格，也是因人而异的。如果在闭合解的范畴以外，再加上一些特别函数(比如Bessel函数、Gamma函数、无限级数、连续分数)等，就是所谓的解析解(Analytical solution)。解析解中不包含极限和积分运算。

(二) 3个概念的整体辨析

解析解

解析解(Analytical solution) 就是根据严格的公式推导，给出任意的自变量就可以求出其因变量，也就是问题的解，然后可以利用这些公式计算相应的问题。所谓的解析解是一种包含分式、三角函数、指数、对数甚至无限级数等基本函数的解的形式。用来求得解析解的方法称为解析法(Analytical techniques)，解析法即是常见的微积分技巧，例如分离变量法等。
封闭解

解析解是一个封闭形式(Closed-form) 的函数，因此对任一自变量，我们皆可将其带入解析函数求得正确的因变量。因此，解析解也被称为封闭解(Closed-form solution)。
数值解

数值解(Numerical solution) 是采用某种计算方法，如有限元法，数值逼近法，插值法等得到的解。别人只能利用数值计算的结果，而不能随意给出自变量并求出计算值。

当无法藉由微积分技巧求得解析解时，这时便只能利用数值分析的方式来求得其数值解了。在数值分析的过程中，首先会将原方程加以简化，以利于后来的数值分析。例如，会先将微分符号改为差分(微分的离散形式)符号等，然后再用传统的代数方法将原方程改写成另一种方便求解的形式。这时的求解步骤就是将一自变量带入，求得因变量的近似解，因此利用此方法所求得的因变量为一个个离散的数值，不像解析解为一连续的分布，而且因为经过上述简化的操作，其正确性也不如解析法可靠。

简而言之：（1）解析解就是给出解的具体函数形式，从解的表达式中就可以算出任何对应值；（2）数值解就是用数值方法求出近似解，给出一系列对应的自变量和解。

0.3 启发式算法的概念

启发式算法(heuristic algorithm)可以这样定义：一个基于直观或经验构造的算法，在可接受的花费（指计算时间和空间）下给出待解决组合优化(NP-Hard)问题每一个实例的一个可行解，该可行解与最优解的偏离程度一般不能被预计。现阶段，启发式算法以仿自然体算法为主，主要有蚁群算法、模拟退火法、神经网络等。

在数学优化和计算机科学中，启发式算法用于在经典方法效果不佳时更快地解决问题或找到近似解。这是通过牺牲最优性、完整性、准确性或精度换取速度来实现的。它是一个与最优化算法相对的概念。

在某种程度上，它可以被认为是一条捷径。因此启发式算法给出的答案是具有偶然性的，或者说时好时坏，因为启发式方法仅仅告诉我们该如何去找答案，而不直接告诉我们答案是什么。

某些启发式算法具有强大的基础理论，它们要么是从理论中以自上而下的方式推导出来的，要么是基于实验或现实世界的数据得出的。其他的只是基于现实世界观察或经验的经验法则，甚至连理论都没有，这些启发式算法的难点是建立符合实际问题的一系列启发式规则。

关于启发式算法更多更详细内容请参考：启发式算法(cha138.com)

0.4 最优化涉及的一些基本概念

(一) 线性规划中几个概念辨析

最优化线性规划问题中的常见概念辨析：可行解，最优解，基，基向量，非基向量，基变量，非基变量等等

线性规划里面有很多基本的概念容易弄混

已知标准型为：

$\begin{cases} \min & \boldsymbol z = \boldsymbol{C}^{\mathrm T} \boldsymbol x \\ \text { s.t. } & \boldsymbol{Ax} = \boldsymbol b \\ & \boldsymbol x \geq \boldsymbol 0 \end{cases}$

其中，$\boldsymbol C \in \mathbb R^n, \boldsymbol x \in \mathbb R^n, \boldsymbol b \in \mathbb R^m, \boldsymbol A \in \mathbb R^{m \times n}(m <n)$，且$\boldsymbol A$的行向量线性无关。

概念	释义
可行解	只要满足约束条件，$ \boldsymbol{Ax} = \boldsymbol b , \boldsymbol x \geq \boldsymbol 0$的解$\boldsymbol x$称为线性规划问题的可行解。
最优解	不仅要满足约束条件，而且要使目标函数$\boldsymbol z = \boldsymbol{Cx}$达到最小值的可行解称为最优解。
基	$\boldsymbol A_{m\times n}$是约束系数矩阵，秩为$m$，若$\boldsymbol B_{m \times m}$是$\boldsymbol A$的子阵且可逆(即$\boldsymbol B$是$\boldsymbol A$中$m$个线性无关向量组)，称$\boldsymbol B$为LP的一个基，不失一般性，可设$\boldsymbol B$是$\boldsymbol A$的前$m$列，即$\boldsymbol B = (\boldsymbol a_1, \boldsymbol a_2, \cdots, \boldsymbol a_m)$，其相对应的变量$\boldsymbol x_B = (x_1, x_2, \cdots, x_m)^{\mathrm T}$为基变量，基变量有$m$个。其余变量$\boldsymbol x_N = (x_{m+1}, x_{m+2}, \cdots, x_n)^{\mathrm T}$称为非基变量，非基变量有$n-m$个。
基向量	基$\boldsymbol B$的一列就称为一个基向量。
非基向量	在$\boldsymbol A$中除了基$\boldsymbol B$之外的一列就称为基$\boldsymbol B$的非基向量。
基本解( 基解基础解)	若在约束方程组系数矩阵中找到一个基，令其非基变量为零，再求解该$m$元线性方程组可得到唯一解，该解称之为线性规划的基本解，其形式为$\boldsymbol x = (x_1, x_2, \cdots, x_m, 0, 0, \cdots, 0)^{\mathrm T}$ 。
基可行解	基解可正可负，负则不可行(违背非负性约束条件)，称满足所有约束条件的基解为基可行解。
退化的基可行解	若某个基变量取值为零，则称之为退化的基可行解。

所以需要注意的是基本解不一定是可行解，非负的基解才是可行解。

说了半天没有一个直观的理解，下面一张图展示了什么是基本解，什么是基本可行解。

基本解：各个等式约束直线的交点，外加与坐标轴的交点；
基本可行解：基本解里面在可行域范围的那些基本解，可行域的顶点；
最优解：基本可行解里面使目标函数最大(最小)的基本可行解

(二) 数学中的「平凡」「非平凡」「退化」等概念

首先，平凡解就是显而易见、没有讨论的必要但为了结果的完整性仍要考虑的解。

比如$\boldsymbol{Ax} = \boldsymbol 0$中的一个显而易见的解——$\boldsymbol x = \boldsymbol 0$，即为平凡解。
比如求一个数的因子，正负的1和它本身是这个数最显而易见的解，所以这两个因子就是平凡解，其他的是不平凡的。

其次，「退化」在数学里面有“不可逆”、“不满秩”这样的含义。一个复杂的东西在某种条件下变成一个简单的东西，就叫「退化」。它一般跟“某个量等于0”有关。① 比如行列式等于0可以叫退化情形(不可逆矩阵)；② 比如三大二阶线性偏微分方程分类里面，抛物型方程一般被称为退化情形，因为它对应的二次型是退化的——实际上也就是说对应的系数矩阵是不可逆的；③ 比如在某一点处的概率为1而在其他点的概率都为0的概率分布，它就比其他一般的分布都更简单，因此叫做退化分布； ④ 比如一个一元二次方程，在二次项系数为零时，就退化成一元一次方程；⑤ 比如椭圆在离心率为零时退化成圆；⑥ 比如向量的内积，在向量只有一维的时候就退化成标量乘法。

退化另一释义一般指是某种特殊，极限的情形。比如说欧式几何里两条平行的直线，比如零乘一个多项式变成零，比如一个环里有两个非零的元素乘积是零，比如坐标变换在一点上雅可比行列式为零，比如原本线性无关的基变换后线性相关了……

补充，机器学习中加正交约束避免平凡解怎么理解呢？

机器学习中很多参数的求解是迭代解法，根据奥卡姆剃刀原理，在结果相同的情况下，模型倾向于选择一个最简单的解。比如：$\boldsymbol{Ax} + \boldsymbol{By} = \boldsymbol 0$，待求参数为$\boldsymbol x$和$\boldsymbol y$，不得不承认$\boldsymbol{x=y=0}$是该方程的解(其实就是上面说的平凡解)，但是这并不是我们想要的结果，因此可以通过增加一个约束项，例如：$|\boldsymbol x|^2>ε$来限制解的范围。平凡解的确也是一种解，但是并不是我们所期望的解，为了避免在迭代过程中模型收敛到平凡解，因此通过增加适当的约束项，来约束解的范围。

(三) 小概念之凸包

平面上给定若干点，则它们的凸包被定义为所有能包含所有点的凸多边形的交集。

在一切包含所有点的凸多边形中，凸包的面积是最小的(显然)，周长也是最小的(因为每个包含所有点的凸多边形，总可以通过下面这样的过程，不断缩小周长，最终得到凸包，完整证明可参见此网页)。类似地可以证明，凸包的每个顶点都是原有的点。

凸包的严格定义为：设集合$S ⊂ \mathbb R^n$是由$\mathbb R^n$中的$k$个点所组成的集合，即$S = (\boldsymbol x_1,\boldsymbol x_2, \cdots, \boldsymbol x_k)$。定义$S$的凸包为：

$\mathrm{conv}(S) = \left\{\boldsymbol x = \sum_{i=1}^k \lambda_i \boldsymbol x_i \left| \sum_{i=1}^k \lambda_i=1, \lambda_i \geq 0 \right.\right\}$

详细内容请参考如下链接：

参考链接1：对凸包定义的可视化理解 - Zewbie的文章 - 知乎
参考链接2：凸包 - 维基百科，自由的百科全书 (wikipedia.org)
参考链接3：算法学习笔记(65): 凸包 - Pecco的文章 - 知乎

优化理论的系统笔记内容：最优化理论与方法 - 快乐江湖的专栏 - CSDN：强推！！！

0.5 DOA估计之确定性VS随机性最大似然法

(一) 概述

在空间谱估计中，根据入射信号的模型，最大似然算法基本上分为两类：确定性最大似然(DML)和随机性最大似然(SML)，具体为：

当信号模型是未知的确定模型时导出的最大似然算法称为DML算法；
当入射信号服从高斯随机分布模型时导出的最大似然算法是SML算法；

窄带远场信号的 DOA数学模型为：

$\boldsymbol X(t) = \boldsymbol A(\boldsymbol \theta) \boldsymbol s(t) + \boldsymbol N(t), \quad t=1,2,\cdots,L$

不同的快拍数之间的噪声协方差矩阵为$\boldsymbol O$，也就是不相关，且各阵元接收的噪声是正态分布的，噪声功率为$\sigma^2$。

(二) DML模型

根据以上的假设可知，对于DML由于模型未知，所以观察数据的一阶和二阶矩满足如下条件：

$\begin{aligned} 均值: \quad & \mathbb E[\boldsymbol x(t_i)] = \boldsymbol A(\boldsymbol \theta) \boldsymbol s(t_i) = \bar{\boldsymbol x}(t_i) \\ 协方差: \quad & \mathbb E[(\boldsymbol x(t_i) - \bar{\boldsymbol x}(t_i)) \cdot (\boldsymbol x(t_i) - \bar{\boldsymbol x}(t_i))^{\mathrm H}] = \sigma^2 \boldsymbol I \delta_{ij}\\ 伪协方差: \quad & \mathbb E[(\boldsymbol x(t_i) - \bar{\boldsymbol x}(t_i)) \cdot (\boldsymbol x(t_i) - \bar{\boldsymbol x}(t_i))^{\mathrm T}] = \boldsymbol O \end{aligned}$

其中，伪协方差的概念请参考复随机过程概念补充 - 博客侦探 - 博客园。根据上式可以得到DML的联合概率密度函数，显然有观测矢量$L$次快拍联合(条件)概率密度函数(PDF)：

$\boldsymbol f_{\mathrm{DML}}(\boldsymbol x_1, \cdots, \boldsymbol x_L) = \prod_{i=1}^L \dfrac{1}{\mathrm{det}(\pi \sigma^2 \boldsymbol I)} \exp\left(-\dfrac{1}{\sigma^2} |\boldsymbol x_i - \boldsymbol{As}_i|^2\right)$

(三) SML模型

对于SML，由于其信号的模型为高斯随机分布，所以观察数据是零均值，其二阶矩满足：

$\begin{aligned} 均值: \quad & \mathbb E[\boldsymbol x(t_i)] = \boldsymbol 0 \\ 协方差: \quad & \mathbb E[\boldsymbol x(t_i) \cdot \boldsymbol x^{\mathrm H}(t_i)] = [\boldsymbol A(\boldsymbol \theta) \boldsymbol R_{\boldsymbol s} \boldsymbol A^{\mathrm H}(\boldsymbol \theta) + \sigma^2 \boldsymbol I] \delta_{ij}\\ 伪协方差: \quad & \mathbb E[\boldsymbol x(t_i) \cdot \boldsymbol x^{\mathrm T}(t_i)] = \boldsymbol O \end{aligned}$

所以单次观察数据的似然函数为：

$f_{i}(\boldsymbol x) =\dfrac{1}{\pi^M \mathrm{det}(\boldsymbol R)} \exp\left(\boldsymbol x^{\mathrm H}_i \boldsymbol{R_s}^{-1} \boldsymbol x_i\right)$

于是很容易得出观测矢量$L$次快拍联合(条件)概率密度函数PDF：

$\boldsymbol f_{\mathrm{SML}}(\boldsymbol x_1, \cdots, \boldsymbol x_L) = \prod_{i=1}^L \dfrac{1}{\pi^M \mathrm{det}(\boldsymbol R)} \exp\left(\boldsymbol x^{\mathrm H}_i \boldsymbol{R_s}^{-1} \boldsymbol x_i\right)$

0.6 待定

为什么核范数能凸近似矩阵的秩？为什么核范数是凸的？ - 知乎
 核范数 Nuclear Norm - 人工智能百科 - 超神经 (hyper.ai)

1. 压缩感知基础

压缩感知(Compressed Sensing, CS)是由陶哲轩等人提出的一种用于信息获取的突破性理论。该理论指出：对于稀疏信号或可压缩信号，可采用低于奈奎斯特采样频率的方式对数据采样，降低数据传输量，并能以高概率精确地重建该信号。其原始文献为：D. L. Donoh, Compressed sensing, IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306.

1.1 压缩感知的流程

压缩感知主要分为三步：

1、信号的稀疏化表示，实现信号的压缩；
2、观测矩阵的设计，得到观测数值；
3、信号的重构，得到恢复信号。

1.2 信号稀疏化表示

稀疏信号定义：设一维离散信号$\boldsymbol x$，长度为$N$，可看作为$N$维空间$N\times1$的列向量，若此列向量中含有$K$个不为0元素，且$K \ll N$，则称该信号$\boldsymbol x$是$K$-稀疏信号，具有稀疏性。$K$称为信号$\boldsymbol x$的稀疏度。牢记，信号具有稀疏性（在本域或其他变换域皆可）是可以使用压缩感知的前提！

如果信号稀疏，则信号$\boldsymbol x$可以表示为：

$\boldsymbol x = \sum_{k=0}^{N-1} \boldsymbol \psi_k s_k = \boldsymbol {\Psi s}$

其中，$\boldsymbol s$为稀疏系数，也是我们使用算法恢复出来的重建新号，其尺寸与信号$\boldsymbol x$相同；$\boldsymbol \Psi$是我们所熟知的稀疏矩阵，其尺寸为$N\times N$，这就是信号的稀疏过程。

1.3 观测矩阵设计

观测矩阵的作用主要是使人们可以看到由仪器所获得的观测值$ \boldsymbol y$，其中，你想要看到多少，这就是我们所熟知的采样率了。对观测矩阵的要求是，可从观测值$\boldsymbol y$中高精度的重构出长度为$N$的原始信号$\boldsymbol s$，或者重构出在稀疏矩阵下的等价信号。其具体表达式为：

$\boldsymbol y = \boldsymbol{\Phi x}$

式中，$\boldsymbol \Phi$为观测矩阵(也称测量矩阵)，其维度为$M\times N$，$\boldsymbol x$为原始信号，其维度为$N\times 1$；$\boldsymbol y$为观测信号，其维度为$M\times 1$。从原始信号中可以观察到多少信号，就是由测量矩阵决定的，因此其会涉及到采样率的概念。通常的采样率概念定义为测量矩阵的行数/列数，即$\dfrac{M}{N}$。

将这个式子带入1.2节中的式子，则可以得到大家经常看到的式子了：

$\boldsymbol y = \boldsymbol{\Phi x} = \boldsymbol{\Phi \Psi s} = \boldsymbol{\Theta s}$

其中，$\boldsymbol {\Theta=\Phi \Psi}$为$M\times N$阶矩阵，又称为传感矩阵。

在此步中，最重要的是构造出合适的观测矩阵，使得可通过仪器采集到从稀疏信号中获得的观测值，并在反向求解时，由观测值重构稀疏信号，即构造出有解的$M \times K$线性方程组。牢记，在进行上述过程时，一定要对观测矩阵进行RIP性质分析，说白了就是观测矩阵和稀疏矩阵的相关性很小很小很小！

1.4 信号重构

信号重构是对$\boldsymbol{y =\Phi x = \Phi \Psi s = \Theta s}$式求最优解，是压缩感知理论中的求解问题，如何得到最优解是研究的主要内容，也是最后一个关键的步骤。压缩感知目前的重构算法主要分为两类：贪婪算法与凸优化算法。

贪婪算法主要是选择合适的列向量经过多次的逐步加和以实现信号的逼近，其中匹配追踪算法、正交匹配追踪等算法均属于贪婪算法；
凸优化算法则是将范数的求解置于范数进行线性规划求解，此算法包括基追踪算法、梯度投影算法等。

最后总结一下：对信号进行压缩感知恢复重建时，信号要映射到稀疏域，即信号在一定的变换域中是稀疏的。在进行测量矩阵压缩采样时，测量矩阵需满足RIP条件，或者应满足观测矩阵与稀疏基不相关的条件。在两个条件的基础上，可运用压缩感知算法进行恢复重建。

2 压缩感知的常用稀疏基

2.1 稀疏基的概念

稀疏基是将可以信号进行稀疏表示的一个矩阵。设一维离散信号$\boldsymbol x$，长度为$N$，可看作为$N$维空间$N \times1$的列向量，若此列向量中含有$K$个不为0元素，且$K \ll N$，则称该信号$\boldsymbol x$是$K$-稀疏信号，具有稀疏性。

如果信号稀疏，则信号$\boldsymbol x$可以表示为：

$\boldsymbol x = \boldsymbol{\Psi s}$

其中，$\boldsymbol \Psi$就是稀疏基(稀疏矩阵)。

2.2 一些常用的稀疏基

DFT稀疏基
DCT(离散余弦)稀疏基
DWT稀疏基
0-1随机稀疏基

3 贪婪算法

压缩感知第三步是进行信号的重构，需要用到恢复重构算法。前面的文章提到过，压缩感知的恢复算法主要分为贪婪算法和凸优化算法两种，这里主要介绍贪婪算法的两种基础算法：MP算法和OMP算法

3.1 匹配追踪算法(MP)

算法假定输入信号与字典库中的原子在结构上具有一定的相关性，这种相关性通过信号与原子库中原子的内积表示，即内积越大，表示信号与字典库中的这个原子的相关性越大，因此可以使用这个原子来近似表示这个信号。

当然这种表示会有误差，将表示误差称为信号残差，用原信号减去这个原子，得到残差，再通过计算相关性的方式从字典库中选出一个原子表示这个残差。迭代进行上述步骤，随着迭代次数的增加，信号残差将越来越小，当满足停止条件时终止迭代，得到一组原子，及残差，将这组原子进行线性组合就能重构输入信号。

3.2 正交匹配追踪算法(OMP)

正交匹配追踪算法迭代的基本思想就是每次迭代过程中从全息矩阵$\boldsymbol T$中选出与测量信号$\boldsymbol s$相关度(内积)最大的那一列，然后从$\boldsymbol T$中去掉该列并加入到扩充矩阵$\boldsymbol {\mathrm{Aug}}_t$中，然后利用最小二乘法原理求出使残差$\boldsymbol r_n =\boldsymbol s-\boldsymbol {\mathrm{Aug}}_t*\boldsymbol {\mathrm{Aug}}_y$最小的一个估计$\boldsymbol {\mathrm{Aug}}_y$，然后不断的从$\boldsymbol T$中减去相关列重复以上过程，直到达到迭代次数结束。

4 凸优化算法

除了贪婪算法以外，压缩感知重构算法的另一大类就是凸优化算法，这类方法通过将非凸问题转化为凸问题求解找到信号的逼近。

4.1 基追踪算法（BP）

该算法全称为Basis Pursuit，其提出使用$l_1$范数替代$l_0$范数来解决最优化问题，以便使用线性规划方法来求解。即将求解 $\min \limits_{\boldsymbol \alpha} ||\boldsymbol \alpha||_{l_{0}}$ 的问题转化为求解 $\min\limits_{\boldsymbol \alpha} ||\boldsymbol \alpha||_{l_{1}}$ 的问题，其中$\boldsymbol \alpha$是需要构建出来的原始信号(重构信号)。

参考链接4.1.1：如何理解“L1范数是L0范数的最优凸近似”？ - 牛狮的回答 - 知乎
参考链接4.1.2：Basis Pursuit, LASSO, TV三种压缩感知优化求解算法 - SINK IN的文章 - 知乎
参考链接4.1.3：压缩感知中l1正则化最小二乘问题的算法研究 - 豆丁建筑

4.2 梯度投影法

梯度投影法(gradient projection method)利用梯度的投影技巧求约束非线性规划问题最优解的一种方法，求带线性约束的非线性规划问题更为有效。它是从一个基本可行解开始，由约束条件确定出凸约束集边界上梯度的投影，以便求出下次的搜索方向和步长。每次搜索后，都要进行检验，直到满足精度要求为止。具体原理可参考梯度投影法 - 百度百科。

本文转自：
版权声明：本文为CSDN博主「爱学习的一一」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处
及本声明。
原文链接：https://blog.csdn.net/weixin_42105848/article/details/124280903

参考链接1：形象易懂讲解算法II——压缩感知 - 咚懂咚懂咚的文章 - 知乎：压缩感知的简要概述。
参考链接2：如何理解压缩感知(compressive sensing)？ - 知乎：压缩感知的简要概述。
参考链接3：压缩感知基本理论-飞大圣的博客-CSDN：压缩感知的简要概述。
参考链接4：随笔分类 - 压缩感知compressive sensing-AndyJee的博客-博客园：强推！！！博客园博主AndyJee关于压缩感知的31篇文章。
参考链接5：Bessel序列与框架的若干刻画 - 豆丁网 (docin.com)
参考链接6：最优化理论——线性规划的标准形式 - 不喝牛奶的文章 - 知乎、优化理论——线性规划 - 星空爱好者的文章 - 知乎
参考链接7：压缩感知 Compressed Sensing - 百科全书、科学新闻和研究评论

5. 原子范数最小化(Atomic Norm Minimization)

原子范数最小化算法能在一众算法中登顶的原因也很简单：它既是拥有优良凸优化性质的算法，又没有精度的限制。简单而言，如果说OMP等在有限码本上选取码字的算法为On-grid类型。那么原子范数最小化算法就是在无穷精度的范围内进行搜索，即Gridless类型。

5.1问题建模

考虑通信中常见的DOA问题。对于$M$根接收天线而言，接收数据可表示为：

$\boldsymbol Y = \boldsymbol A(f) \boldsymbol S + \boldsymbol E$

其中，$\boldsymbol Y \in \mathbb C^{M \times L}$，$L$代表接收时隙数(可理解为观测次数)，$\boldsymbol S \in \mathbb C^{K \times L}$代表了这$K$个源在$L$个时隙的发射信号，$\boldsymbol E \in \mathbb C^{M \times L}$为噪声，$\boldsymbol A(f)$为：

$\boldsymbol A(f) = [\boldsymbol a(f_1), \cdots, \boldsymbol a(f_K)] \in \mathbb C^{M \times K}$

其每一列对应第$k$个源(source)的DOA对应的天线响应矢量，式中：

$\boldsymbol a(f_i) = \left[1 \quad \mathrm e^{\mathrm j 2 \pi f_i} \quad \cdots \quad \mathrm e^{\mathrm j 2 \pi(M-1) f_i}\right]^{\mathrm T}$

其中，$f_i = \cos \theta_i / 2$。

Toeplitz矩阵的范德蒙德分解 ：
对于任意的秩满足$r \leq N$的半正定Toeplitz矩阵$\boldsymbol T(\boldsymbol u) \in \mathbb{C}^{N \times N}$，则有如下的$r-$原子范德蒙德分解： $$ \boldsymbol T(\boldsymbol u) = \sum_{k=1}^{r} p_k \boldsymbol a(f_k) \boldsymbol a^{\mathrm H}(f_k) = \boldsymbol A(f) \mathrm{diag}(\boldsymbol p) \boldsymbol A^{\mathrm H}(f) $$ 其中，$\boldsymbol A(f) = [\boldsymbol a(f_1), \boldsymbol a(f_2), \cdots, \boldsymbol a(f_r)]$。当$r < N$时，此分解是唯一的。
此定律十分重要，原子范数的理论推导直接建立在这个分解之上！！！

5.2 DOA估计的一般框架

考虑$L = 1$即一次观测的简单场景，表示为：

$\boldsymbol y = \boldsymbol z + \boldsymbol e = \boldsymbol A(f) \boldsymbol s + \boldsymbol e$

那么一般可以以如下的优化问题为目标：

$\begin{cases} \min & \mathcal{M}(\boldsymbol{z}) \\ \text { s.t. } & ||\boldsymbol{y - z}||_2 \leq \varepsilon \end{cases}$

其中，$\mathcal{M}(\boldsymbol z)$是选择的一种度量metric。进一步地，通过引入惩罚系数$\lambda$，可将其改写为：

$\min_{\boldsymbol z} \lambda \mathcal{M}(\boldsymbol z) + \dfrac{1}{2} ||\boldsymbol{y - z}||_2$

下面要重点介绍的，就是以原子范数$||\boldsymbol z||_{\mathcal A}$作为$\mathcal{M}(\boldsymbol z)$的方法。

5.3 原子范数的概念

(一) $l_0$原子范数

首先，原子集的定义如下：

$\mathcal A = \left\{\boldsymbol a(f, \phi) = \boldsymbol a(f) \phi \mid f \in \mathbb T, \phi \in \mathbb C, |\phi| = 1 \right\}$

这个集合可以理解为是类似于OMP方法中的字典，但它是无限精度的，因为$f$可以是任意实数。$\phi$是允许了一个初始相位的不同。而根据$\boldsymbol z = \boldsymbol A(f) \boldsymbol s$，显然$\boldsymbol z$是该原子集中$k$个原子的线性组合。而$\ell_0$-原子范数，就是指能组成$\boldsymbol z$的最少所需原子数，即：

$\begin{aligned} ||\boldsymbol z||_{\mathcal A, 0} =& \inf_{c_k, f_k, \phi_k} \left\{ \mathcal K \left| \boldsymbol z = \sum_{k=1}^{\mathcal K} c_k \boldsymbol a(f_k, \phi_k), \quad c_k >0, f_k \in \mathbb T, |\phi_k| = 1 \right.\right\} \\ = & \inf_{s_k, f_k} \left\{ \mathcal K \left| \boldsymbol z = \sum_{k=1}^{\mathcal K} s_k \boldsymbol a(f_k), \quad f_k \in \mathbb T \right. \right\} \end{aligned}$

因为我们的目的就是为了恢复出$\boldsymbol A(f)$，而$\boldsymbol z$可以写成无数种$\mathcal{A}$ 中原子线性组合的形式。但只有对应所用原子数最少即对应最小$\ell_0$-原子范数时，此时组成$\boldsymbol z$的原子才恰好对应待恢复的$\boldsymbol A(f)$。因此，下面的任务就是最小化$\boldsymbol z$的$\ell_0$-原子范数$||\boldsymbol z||_{\mathcal A, 0}$。

(二) $\ell_0$-原子范数与范德蒙德分解

虽然模型建立出来了，但是容易看到，如何对$||\boldsymbol z||_{\mathcal A, 0}$进行最小化，可以说是完全摸不着头脑。这似乎比$\ell_0$范数最小化更为抽象。此时就是见证数学魅力的时刻，考虑如下形式的矩阵：

$\left[\begin{array}{cc} x & \boldsymbol z^{\mathrm H} \\ \boldsymbol z & \boldsymbol T(\boldsymbol u) \end{array}\right] \succeq 0$

$\boldsymbol T(\boldsymbol u)$的定义如之前，是一个由$\boldsymbol u$得到的Teoplitz矩阵。$x$则是一个待优化的变量。这个约束隐含了如下的结论：

① $\boldsymbol T(\boldsymbol u) \succeq 0$，否则无法保证对于任何向量$\boldsymbol y$，都有$\boldsymbol y^{\mathrm H} \left[\begin{array}{cc} x & \boldsymbol z^{\mathrm H} \\ \boldsymbol z & \boldsymbol T(\boldsymbol u) \end{array}\right]\boldsymbol y \geq 0$
② $\boldsymbol z$ 一定位于$\boldsymbol T(\boldsymbol u)$的列空间中。这个证明也在后面的证明章节中给出了。

由于$\boldsymbol T(u) \succeq 0$，因此$\boldsymbol T(u)$是一个半正定矩阵，即存在范德蒙德分解：

$\boldsymbol T(u) = \sum_{k=1}^r p_k \boldsymbol a(f_k) \boldsymbol a^{\mathrm H}(f_k) = \boldsymbol A(f) \mathrm{diag}(\boldsymbol p) \boldsymbol A^{\mathrm H}(f)$

其中，$r = \mathrm{rank}(\boldsymbol T(u))$。由于$\boldsymbol z$一定位于$\boldsymbol T(u)$的列空间中，那么$\boldsymbol z$必能写为$\boldsymbol a(f_k)$的线性组合！这一点至关重要，因为这引出了如下的结论：

最小化$\ell_0$-原子范数等价于求解如下问题：

$\begin{cases} \min\limits_{x, u} & r = \mathrm{rank}(\boldsymbol T(u)) \\ \text { s.t. } & \left[\begin{array}{cc} x & \boldsymbol z^{\mathrm H} \\ \boldsymbol z & \boldsymbol T(u) \end{array}\right] \succeq 0 \end{cases}$

由于$\boldsymbol z$是$\boldsymbol T(u)$范德蒙德分解所得的$r$个原子的线性组合。当找到秩最小的$r$时，也就找到了组成$\boldsymbol z$所需的最少原子数。也就对应$\boldsymbol z$的$\ell_0$-原子范数最小化。而对于这个转化后的问题，如果$u$被解出，那么$\boldsymbol T(u)$也能得到，那么对$\boldsymbol T(u)$进行范德蒙德分解，也就获得了$\boldsymbol A(f)$。但是美中不足的是，该目标函数并非凸函数，无法轻易求解。

(三) $\ell_1$-原子范数

将$\ell_0$-原子范数进行凸松弛，得到的就是$\ell_1$-原子范数。其定义如下：

$||\boldsymbol z||_{\mathcal A, 1} = \inf_{f_k, s_k} \left\{ {\sum_k ∣s_k∣ \left| \boldsymbol z = \sum_k \boldsymbol a(f_k) s_k , \quad f_k \in \mathbb T \right.} \right\}$

与$\ell_0$-原子范数的定义进行比较，发现这和将传统的$\ell_0$-范数松弛为$\ell_1$-范数如出一辙。然而如何最小化$||\boldsymbol z||_{\mathcal A, 1}$看上去也非常困难。此时再度利用范德蒙德分解，有如下精彩的结论。

最小化$\ell_1$-原子范数等价于求解如下问题：

$\begin{cases} \min\limits_{x, u} & \dfrac{1}{2}x + \dfrac{1}{2}u \\ \text { s.t. } & \left[\begin{array}{cc} x & \boldsymbol z^{\mathrm H} \\ \boldsymbol z & \boldsymbol T(u) \end{array}\right] \succeq 0 \end{cases}$

注意到，这是一个凸问题。首先目标函数显然是变量的仿射函数，因此为凸(既凸且凹)。限制条件也可以写为变量的仿射函数形式。因此也满足凸问题的限制条件($f(x) \leq 0$，$f(x)$为凸函数)。($\boldsymbol X \succeq 0$可以等价为$\boldsymbol y^{\mathrm H} \boldsymbol X \boldsymbol y \geq 0, \quad \forall \boldsymbol y$，而对于每个$\boldsymbol y$，都是关于$\boldsymbol X$的仿射变换）

$\boldsymbol T(u)$的对角元素都是$u$，那么事实上$u$就是$\mathrm{tr}(\boldsymbol T(u))$！而后者又被称为核函数，也是$\mathrm{rank}(\boldsymbol T(u))$的经典凸松弛。这从另一个角度解释了$\ell_1$原子范数最小化是$\ell_0$原子范数最小化的凸松弛。

至此，压缩感知问题被转化为了一个可以由CVX进行直接求解的凸问题！还剩的最后一块拼图：即为何$\ell_1$-原子范数可以等效为这个凸问题，仍照例，放在下面的证明章节中。

(四) 多维原子范数

将原子范数拓展到多维是十分必要的，因为通信中DOA估计大多是多次观测。然而其结论大体相似。此时，变量变为二维矩阵$\boldsymbol Z$，而其$\ell_0$-原子范数被定义为：

$\begin{aligned} ||\boldsymbol Z||_{\mathcal A, 0} =& \inf_{c_k, f_k, \boldsymbol \phi_k} \left\{ \mathcal K \left| \boldsymbol Z = \sum_{k=1}^{\mathcal K} c_k \boldsymbol a(f_k, \boldsymbol \phi_k), \quad c_k >0, f_k \in \mathbb T, ||\boldsymbol \phi_k||_2 = 1 \right. \right\} \\ = & \inf_{\boldsymbol s_k, f_k} \left\{ \mathcal K \left| \boldsymbol Z = \sum_{k=1}^{\mathcal K} \boldsymbol s_k \boldsymbol a(f_k), \quad f_k \in \mathbb T\right.\right\} \end{aligned}$

其中，$\boldsymbol a(f_k, \boldsymbol \phi_k)$为(注意加粗)：

$\boldsymbol a(f_k, \boldsymbol \phi_k) = \boldsymbol a(f_k) \boldsymbol \phi_k, \quad f_k \in \mathbb T, \boldsymbol \phi_k \in \mathbb C^{1 \times L}, || \boldsymbol \phi_k||_2 = 1$

注意到，这和向量$\boldsymbol z$的原子范数的最大区别在于标量$\phi_k$变为了行向量$\boldsymbol \phi_k$。类似的，对其凸松弛为$\ell_1$原子范数，定义为：

$\begin{aligned} ||\boldsymbol Z||_{\mathcal A, 1} =& \inf \left\{ t>0 \left| \boldsymbol Z \in t\mathrm{conv}(\mathcal A)\right. \right\} \\ =& \inf_{c_k, f_k, \boldsymbol\phi_k} \left\{ \sum_k c_k \left| \boldsymbol Z = \sum_{k=1}^{\mathcal K} c_k \boldsymbol a(f_k, \boldsymbol \phi_k), \quad c_k >0, f_k \in \mathbb T, ||\boldsymbol \phi_k||_2 = 1 \right. \right\} \\ = & \inf_{\boldsymbol s_k, f_k} \left\{ \sum_k ||\boldsymbol s_k||_2 \left| \boldsymbol Z = \sum_{k=1}^{\mathcal K} \boldsymbol s_k \boldsymbol a(f_k), \quad f_k \in \mathbb T\right.\right\} \end{aligned}$

其中，$\mathrm{conv}(\mathcal A)$表示原子集合的凸包。而对其的最小化，仍可以等价为如下SDP问题：

$\begin{cases} \min\limits_{\boldsymbol X, u} & \dfrac{1}{2\sqrt{N}} \left[\mathrm{tr}({\boldsymbol X}) + \mathrm{tr}({\boldsymbol T}(u)) \right] \\ \text { s.t. } & \left[\begin{array}{cc} \boldsymbol X & \boldsymbol Z^{\mathrm H} \\ \boldsymbol Z & \boldsymbol T(u) \end{array}\right] \succeq 0 \end{cases}$

证明见5.4。

5.4 相关证明

(一) 范德蒙德分解的证明

() $\boldsymbol z$存在于列空间的证明

当$\boldsymbol T(\boldsymbol u)$满秩时，直接成立。当$\boldsymbol T(\boldsymbol u)$不满秩时，此时，可以看到：

$\left[\begin{array}{cc} x & \boldsymbol z^{\mathrm H} \\ \boldsymbol z & \boldsymbol T(\boldsymbol u) \end{array}\right] \succeq 0 \Longleftrightarrow [ t \quad \boldsymbol y]^{\mathrm H} \left[\begin{array}{cc} x & \boldsymbol z^{\mathrm H} \\ \boldsymbol z & \boldsymbol T(\boldsymbol u) \end{array}\right] \left[\begin{array}{c} \boldsymbol t \\ \boldsymbol y \end{array}\right] \geq 0, \forall t,\boldsymbol y \\ \Downarrow\\ t^*x t + \boldsymbol y^{\mathrm H} \boldsymbol{z} t + t^*\boldsymbol z^{\mathrm H}\boldsymbol y + \boldsymbol y^{\mathrm H} \boldsymbol T(\boldsymbol u) \boldsymbol y \geq 0$

此时，取$\boldsymbol y$位于$\boldsymbol T(\boldsymbol u)$的零空间中，则上式进一步变为：

$t^*x t + \boldsymbol y^{\mathrm H} \boldsymbol{z} t + t^*\boldsymbol z^{\mathrm H}\boldsymbol y \geq 0$

若$\boldsymbol z$不位于$\boldsymbol T(\boldsymbol u)$的列空间中，即$a = \boldsymbol y^{\mathrm H} \boldsymbol z \neq 0$，则：

$t^*x t + a t + t^*a^* \geq 0$

此时，取$t = -\dfrac{a^*}{x}$，则有：

$\dfrac{aa^*}{x^*} - \dfrac{aa^*}{x^*} - \dfrac{aa^*}{x^*} \geq 0$

由于$x$必为非负实数(取$t = 1, y = 0$)，所以上式显然不成立，矛盾。得证。

() 原子范数最小化的等价性证明

记等价凸问题的目标函数为：

$F=\dfrac{1}{2} x+\dfrac{1}{2} u_{1}$

先证明：$F \leq ||\boldsymbol z||_{\mathcal A}$。

取$\boldsymbol z = \sum\limits_{k}c_k \boldsymbol a(f_k, \phi_k) = \sum\limits_{k}s_k \boldsymbol a(f_k)$，令$\boldsymbol T(\boldsymbol u) = \sum\limits_{k} c_k \boldsymbol a(f_k) \boldsymbol a^{\mathrm H}(f_k)$、$x = \sum\limits_k c_k$，则有：

$\left[\begin{array}{cc} x & \boldsymbol{z}^{\mathrm{H}} \\ \boldsymbol{z} & \boldsymbol{T} \end{array}\right]=\sum_{k} c_{k}\left[\begin{array}{c} \phi_{k}^* \\ \boldsymbol{a}\left(f_{k}\right) \end{array}\right]\left[\begin{array}{c} \phi_{k}^* \\ \boldsymbol{a}\left(f_{k}\right) \end{array}\right]^{\mathrm{H}} \succeq 0$

因此$\boldsymbol T$和$x$是一组可行解。那么其必然不小于最优解。而这组可行解对应的目标函数为：

$\dfrac{1}{2} x+\dfrac{1}{2} u_{1} = \sum_k c_k$

因此$F \leq ||\boldsymbol z||_{\mathcal A}$。

再证明：$F \geq ||\boldsymbol z||_{\mathcal A}$

设凸问题的最优解为$(\hat x , \hat u)$，对$\boldsymbol T (\hat u)$进行范德蒙德分解，得到参数 $(\hat r, \hat p_k, \hat f_k)$。根据之前的分析，$\boldsymbol z$必然在$\boldsymbol T(\hat u)$的列空间中，即：

$\boldsymbol z = \sum_{k=1}^{\hat r} \hat c_k \boldsymbol a(\hat f_k, \hat\phi_k) = \sum_{k=1}^{\hat r} \hat s_k \boldsymbol a(\hat f_k)$

进一步地，根据Schur补条件：

$\left[\begin{array}{cc} x & \boldsymbol z^{\mathrm H} \\ \boldsymbol z & \boldsymbol T(u) \end{array}\right] \succeq 0 \quad \Longleftrightarrow \quad \hat x \geq \boldsymbol z^{\mathrm H}[\boldsymbol T(\boldsymbol{\hat u})]^{\dagger} \boldsymbol z = \sum_{k=1}^{\hat r}\dfrac{\hat c^2_k}{\hat p_k}$

而$\hat{u}_{1}=\sum\limits_{k=1}^{\hat{r}} \hat{p}_{k}$，因此

$\begin{aligned} F & =\dfrac{1}{2} \hat{x}+\dfrac{1}{2} \hat{u}_1 \\ & \geq \dfrac{1}{2} \sum_{k} \dfrac{\hat{c}_{k}^2}{\hat{p}_{k}}+\dfrac{1}{2} \sum_k \hat{p}_{k} \\ & \geq \sum_{k} \hat{c}_{k} \\ & \geq\|\boldsymbol{z}\|_{\mathcal{A}} . \end{aligned}$

第一个不等号来自于Schur补条件。第二个不等号来自于$\dfrac{1}{x}+x$的最大化问题。第三个不等式来自于$|| \boldsymbol z||_{\mathcal{A}}$的定义，即所有线性分解中，$\sum\limits_k c_k$的最小值。
至此，得证：$F =|z|_{\mathcal{A}}$。

()多维原子范数最小化的等价性证明

6. 稀疏贝叶斯学习(SBL)

稀疏贝叶斯学习 1 - 本帅哥屏蔽了凡人 - CSDN

压缩感知学习总结及Matlab代码实现 - WinrenJoe - CSDN

压缩感知、一阶范数、Lasso - CSDN

Basis Pursuit, LASSO, TV - SINK IN的文章 - 知乎