原子范数初探：以到达角估计为例

到达方向（Direction-of-arrival, DOA）估计是指从形成传感器阵列的多个接收天线的输出中检索若干电磁波/源的方向信息的过程。DOA估计是阵列信号处理中的一个主要问题，在雷达、声纳、无线通信中有着广泛的应用。

基本数学模型

考虑$K$个窄带的远场源信号$s_k,k=1,\cdots ,K$，多个接收阵列从不同方向 { θ k } k = 1 K {\{\theta_k\}}_{k=1}^K {θk​}k=1K​接收到信号。而不同阵列之间的信号延时可以用简单的相移(phase shift)来表征，这也就引出来下面的信号模型：
$$\mathbf{y}(t)=\sum _{k=1}^K\mathbf{a}({\theta }_k)s_k(t)+\mathbf{e}(t)=\mathbf{A}(\mathbf{\theta })\mathbf{s}(t)+\mathbf{e}(t),t=1,\cdots ,L\text{\hspace{1em}(1)}$$

其中$t$表示第$t$个时隙(snapshot)，$L$表示总的时隙数量(the number of snapshots)。$\mathbf{y}(t)\in {\mathbb{C}}^{M\times 1}$，$\mathbf{s}(t)\in {\mathbb{C}}^{K\times 1}$，$\mathbf{e}(t)\in {\mathbb{C}}^{M\times 1}$分别表示阵列在第$t$个时隙的输出信号，源信号向量，以及测量噪声。$M$是阵列的数量。$\mathbf{a}({\theta }_k)$被称作第$k$个信号源的导向矢量(steering vector)，这$K$个导向矢量构成了阵列流形矩阵(array manifold matrix) $\mathbf{A}(\mathbf{\theta })=\left[\mathbf{a}({\theta }_1),\cdots ,\mathbf{a}({\theta }_K)\right]\in {\mathbb{C}}^{M\times K}$。我们可以把式(1)写为更紧凑的矩阵形式：
$$\mathbf{Y}=\mathbf{A}(\mathbf{\theta })\mathbf{S}+\mathbf{E}\text{\hspace{1em}(2)}$$

给定接收信号$\mathbf{Y}\in {\mathbb{C}}^{M\times L}$，以及映射$\theta \to \mathbf{a}(\theta )$，我们的目标是估计出参数 { θ k } k = 1 K {\{\theta_k\}}_{k=1}^K {θk​}k=1K​（即DOAs）。另外，在实际模型中，我们一般没有$K$的先验知识，且我们一般假设$K\le M$。

阵列的几何结构

我们考虑一个简单的线性阵列(linear array)，如下图所示

$\mathbf{a}(\theta )\in {\mathbb{C}}^{M\times 1}$表示为：
$$a_m(\theta )=e^{i\pi r_m\cos \theta },m=1,2,\cdots ,M\text{\hspace{1em}(3)}$$

我们可以把变量$\theta$替换为$f=\frac{1}{2}\cos \theta$，并且定义(without ambiguity)$\mathbf{a}(f)=\mathbf{a}\left(\arccos (2f)\right)=\mathbf{a}(\theta )$，那么映射$\mathbf{a}(f)$为：
$$a_m(f)=e^{i2\pi r_{m}f},m=1,2,\cdots ,M\text{\hspace{1em}(4)}$$

从上式可以看出，空域上的估计问题变成了一个时频估计（temporal frequency estimation）问题（线谱估计问题，linear spectral estimation）。
对时频域的理解：

1. 时间：注意到$r_m$其实想表征第$m-1$根天线距离第0根天线之间的距离，因为距离与时间的正比关系，所以$r_m$对应到时间(temporal)；
2. 频率：注意到$f$的定义域$f\in (-\frac{1}{2},\frac{1}{2}]$，$2\pi f\in (-\pi ,\pi ]$，类比角频率。我们定义$\mathbb{T}=(-\frac{1}{2},\frac{1}{2}]$（之后会用到）

如果我们考虑等间距的线性阵列，那么该阵列即为ULA( Uniform Linear Array )，并且仅考虑最简单的单位时间，即两个相邻天线的时间差归一，对应到$r_m-r_{m-1}=1$，这时
$$\mathbf{a}(f)=[1,e^{i2\pi f},\cdots ,e^{i2\pi (M-1)f}{]}^T\in {\mathbb{C}}^{M\times 1}\text{\hspace{1em}(5)}$$

如果通过仅保留一部分天线，从ULA获得线性阵列，则称为稀疏线性阵列（SLA，Sparse Linear Array）。

ULA模型和SLA模型

ULA模型：$M$根接收天线，则阵列接收信号模型为：
$$\mathbf{Y}=\mathbf{A}(\mathbf{f})\mathbf{S}+\mathbf{E}\text{\hspace{1em}(6)}$$

其中$f_k\in \mathbb{T}$，$k=1,\cdots ,K$是我们感兴趣的频率成分，它与DOAs有着one-to-one的映射关系。

SLA模型：SLA有$M$根接收天线，若这$M$根天线是从$N$个ULA阵列中选取的，记索引为 Ω = { Ω 1 , ⋯   , Ω M } \Omega=\{ \Omega_1, \cdots, \Omega_M \} Ω={Ω1​,⋯,ΩM​}，其中$N\ge M$，且$1\le {\Omega }_1<\cdots <{\Omega }_M\le N$，那么SLA的数据模型为
$${\mathbf{Y}}_{\Omega }={\mathbf{A}}_{\Omega }(f)\mathbf{S}+{\mathbf{E}}_{\Omega }\text{\hspace{1em}(7)}$$

单时隙简单场景(The Single Snapshot Case)

确定性方法的一般框架

SLA中，单时隙场景的数据模型为：
$${\mathbf{y}}_{\Omega }={\mathbf{z}}_{\Omega }+{\mathbf{e}}_{\Omega },\mathbf{z}=\mathbf{A}(\mathbf{f})\mathbf{s}\text{\hspace{1em}(8)}$$

ULA是 Ω = { 1 , ⋯   , N } \Omega=\{1,\cdots,N\} Ω={1,⋯,N}的特殊情况。对于确定性的稀疏方法，一般的，我们需要解决下面的约束优化问题：
$$\underset{\mathbf{z}}{\min }\mathcal{M}(\mathbf{z}),\text{subject to }{\Vert {\mathbf{z}}_{\Omega }-{\mathbf{y}}_{\Omega }\Vert }_2\le \eta \text{\hspace{1em}(9)}$$

其中噪声是有界的：$\Vert {\mathbf{e}}_{\Omega }{\Vert }_2\le \eta$。$\mathcal{M}(\mathbf{z})$表征一种稀疏的metric：最小化$\mathcal{M}(\mathbf{z})$意味着相应的$\mathbf{a}(f)$的数量（the number of components/atoms $\mathbf{a}(f)$）降低。这些原子(atoms)可以直接给出频率的估计结果。我们可以考虑下面的正则优化问题（regularized optimization problem）:
$$\underset{\mathbf{z}}{\min }\lambda \mathcal{M}(z)+\frac{1}{2}{\Vert {\mathbf{z}}_{\Omega }-{\mathbf{y}}_{\Omega }\Vert }_2^2\text{\hspace{1em}(10)}$$

令$\eta \to 0$，可以得到退化的无噪声问题：
$$\underset{\mathbf{z}}{\min }\mathcal{M}(\mathbf{z}),\text{subject to }{\mathbf{z}}_{\Omega }={\mathbf{y}}_{\Omega }\text{\hspace{1em}(11)}$$

之后的内容，我们主要考虑metric $\mathcal{M}(\mathbf{z})$的选择。

注：这里的一个原子(atoms)指的是一个频率成分，即对应到一个$\mathbf{a}(f)$

原子${\ell }_0$范数(Atomic ${\ell }_0$ norm)

我们定义原子集合(the set of atoms)：
A = { a ( f , ϕ ) = a ( f ) ϕ : f ∈ T , ϕ ∈ C , ∣ ϕ ∣ = 1 } (12) \mathcal A = \{ \boldsymbol a(f, \phi) = \boldsymbol a(f) \phi : f \in \mathbb T, \phi \in \mathbb C, |\phi|=1\} \tag{12} A={a(f,ϕ)=a(f)ϕ:f∈T,ϕ∈C,∣ϕ∣=1}(12)

真实信号$\mathbf{z}$是原子集合$\mathcal{A}$中$K$个原子的线性组合，原子${\ell }_0$范数(the atomic ${\ell }_0$ (pseudo-) norm)，我们用$\Vert \mathbf{z}{\Vert }_{\mathcal{A},0}$来表示：is defined as the minimum numer of atoms in $\mathcal{A}$ that can synthesize $\mathbf{z}$：
∥ z ∥ A , 0 = inf ⁡ c k , f k , ϕ k { K : z = ∑ k = 1 K a ( f k , ϕ k ) c k ,    f k ∈ T , ∣ ϕ k ∣ = 1 , c k > 0 } = inf ⁡ f k , s k { K : z = ∑ k = 1 K a ( f k ) s k , f k ∈ T } (13)

\tag{13} ∥z∥A,0​​=ck​,fk​,ϕk​inf​{K:z=k=1∑K​a(fk​,ϕk​)ck​,  fk​∈T,∣ϕk​∣=1,ck​>0}=fk​,sk​inf​{K:z=k=1∑K​a(fk​)sk​,fk​∈T}​(13)

为了解决该问题，我们需要引入Toepitz协方差矩阵的范德蒙德分解。更具体一些，令$\mathbf{T}(u)$为Toeplitz矩阵，并且满足下述条件：
[ x z H z T ( u ) ] ≽ 0 (14) \left[

\right] \succcurlyeq \boldsymbol{0} \tag{14} [xz​zHT(u)​]≽0(14)
其中$\succcurlyeq 0$表示半正定矩阵。$x$是一个待优化的自由变量，如果式(14)要成立，那么$\mathbf{T}(\mathbf{u})$必然是一个半正定的Toeplitz矩阵，因此可以进行$\text{rank}(\mathbf{T}(\mathbf{u}))$-atomic Vandermonde deomposition。另外，$\mathbf{z}\in \text{span}(\mathbf{T}(\mathbf{u}))$，因此原子${\ell }_0$范数与$\mathbf{T}(\mathbf{u})$的秩密切相关，我们有下面的定理

式(13)定义的$\Vert \mathbf{z}{\Vert }_{\mathcal{A},0}$，等于下面优化问题的最优值：
min ⁡ x , u rank ( T ( u ) ) subject to [ x z H z T ( u ) ] ≽ 0 (15) \min_{\boldsymbol x, \boldsymbol u} \text{rank} \left ( \boldsymbol T(\boldsymbol u) \right) \\ \text{subject to} \left[

xzHzT(u)" role="presentation">xzzHT(u)$$\begin{array}{cc}x& {\mathit{z}}^H\\ \mathit{z}& \mathit{T}\left(\mathit{u}\right)\end{array}$$

\right] \succcurlyeq \boldsymbol{0} \tag{15} x,umin​rank(T(u))subject to[xz​zHT(u)​]≽0(15)

原子范数(Atomic Norm)

比原子${\ell }_0$范数更为实用的metric是原子范数最小化(Atomic Norm Minimization, ANM)，我们正式地引入原子范数，原子范数的实质为式(12)集合中$\mathcal{A}$的凸包(convex hull)：
∥ z ∥ A = inf ⁡ { t > 0 : z ∈ t ⋅ c o n v ( A ) } = inf ⁡ c k , f k , ϕ k { ∑ k c k :    z = ∑ k a ( f k , ϕ k ) c k ,   f k ∈ T ,   ∣ ϕ k ∣ = 1 ,   c k > 0 } = inf ⁡ f k , s k { ∑ k ∣ s k ∣ : z = ∑ k a ( f k ) s k ,    f k ∈ T } (16)

\tag{16} ∥z∥A​​=inf{t>0:z∈t⋅conv(A)}=ck​,fk​,ϕk​inf​{k∑​ck​:  z=k∑​a(fk​,ϕk​)ck​, fk​∈T, ∣ϕk​∣=1, ck​>0}=fk​,sk​inf​{k∑​∣sk​∣:z=k∑​a(fk​)sk​,  fk​∈T}​(16)

根据定义，我们不难发现，原子范数其实对应到有限维向量的${\ell }_1$范数。为了解决原子范数问题，学界和业界提供了一种高效的SDP优化方法。

定理：式(16)中定义的$\Vert \mathbf{z}{\Vert }_{\mathcal{A}}$，等于下面SDP问题的最优值：
min ⁡ x u 1 2 x + 1 2 u 1 subject to  [ x z H z T ( u ) ] ≽ 0 (17)

minxu12x+12u1subject to [xzHzT(u)]≽0" role="presentation">minxu12x+12u1subject to [xzzHT(u)]≽0$$\begin{array}{rl}& \underset{x\mathit{u}}{min}\frac{1}{2}x+\frac{1}{2}{u}_1\\ & \text{subject to }\left[\begin{array}{cc}x& {\mathit{z}}^H\\ \mathit{z}& \mathit{T}\left(\mathit{u}\right)\end{array}\right]\succcurlyeq \mathbf{0}\end{array}$$

\tag{17} ​xumin​21​x+21​u1​subject to [xz​zHT(u)​]≽0​(17)

证明：令$F$为式上述优化问题的最优值，我们需要证明$F=\Vert \mathbf{z}{\Vert }_{\mathcal{A}}$
(1)我们首先证明$F\le \Vert \mathbf{z}{\Vert }_{\mathcal{A}}$
令$\mathbf{z}={\sum }_k{c}_k\mathbf{a}(f_k,{\varphi }_k)={\sum }_k\mathbf{a}(f_k)s_k$表征$\mathbf{z}$的原子分解形式。然后令$\mathbf{u}$满足$\mathbf{T}(\mathbf{u})={\sum }_k{c}_k\mathbf{a}(f_k){\mathbf{a}}^H(f_k)$，且$x={\sum }_k{c}_k$，由此我们可以得到
[ x z H z T ( u ) ] = ∑ k c k [ ϕ k ∗ a ( f k ) ] [ ϕ k ∗ a ( f k ) ] H ≽ 0 (18) \left[

\right] =\sum_k{c_k\left[ \right]}\left[ \right] ^H \succcurlyeq \boldsymbol{0} \tag{18} [xz​zHT(u)​]=k∑​ck​[ϕk∗​a(fk​)​][ϕk∗​a(fk​)​]H≽0(18)

因此，由(18)式构建的$x$和$\mathbf{u}$构成了问题(17)的一个可行解，并且目标函数满足：
$$\frac{1}{2}x+\frac{1}{2}{u}_1=\sum _k{c}_k\text{\hspace{1em}(19)}$$

由此，我们可以得到：$F\le {\sum }_k{c}_k$，因为该不等式对$\mathbf{z}$的任意原子分解都成立，因此$F\le {\sum }_k{c}_k$必然成立。

(2)我们证明$F\ge \Vert \mathbf{z}{\Vert }_{\mathcal{A}}$
假设$(\hat{x},\hat{\mathbf{u}})$是问题(17)的一个最优解。因为$\mathbf{T}(\hat{\mathbf{u}})\succcurlyeq 0$，利用半正定Toeplitz矩阵的范德蒙德分解中的定理，我们得到$\mathbf{T}(\hat{\mathbf{u}})$可以进行范德蒙德分解，将分解结果记为$(\hat{r},{\hat{p}}_k,{\hat{f}}_k)$。另外，因为 [ x ^ z H z T ( u ^ ) ] ≽ 0 \left[

\right] \succcurlyeq \boldsymbol{0} [x^z​zHT(u^)​]≽0，可以得到$\mathbf{z}\in \text{span}(\mathbf{T}(\hat{\mathbf{u}}))$，又因为$\text{span}(T(\hat{\mathbf{u}}))=\text{span}\left(\mathbf{a}(f_1),\cdots ,\mathbf{a}(f_{\hat{r}})\right)$，我们知道
$$\mathbf{z}=\sum _{k=1}^{\hat{r}}{\hat{c}}_k\mathbf{a}({\hat{f}}_k,{\hat{\varphi }}_k)=\sum _{k=1}^{\hat{r}}\mathbf{a}({\hat{f}}_k){\hat{s}}_k\text{\hspace{1em}(20)}$$

又因为 [ x ^ z H z T ( u ^ ) ] ≽ 0 \left[

\right] \succcurlyeq \boldsymbol{0} [x^z​zHT(u^)​]≽0，

x ^ ≥ z H [ T ( u ^ ) ] † z = ∑ k = 1 r ^ c ^ k 2 p ^ k u ^ 1 = ∑ k = 1 r ^ p ^ k (21)

\tag{21} x^u^1​​≥zH[T(u^)]†z=k=1∑r^​p^​k​c^k2​​=k=1∑r^​p^​k​​(21)

由此可以得到
F = 1 2 x ^ + 1 2 u ^ 1 ≥ 1 2 ∑ k c ^ k 2 p ^ k + 1 2 ∑ k p ^ k ≥ ∑ k c k ≥ ∥ z ∥ A (22)

\tag{22} F​=21​x^+21​u^1​≥21​k∑​p^​k​c^k2​​+21​k∑​p^​k​≥k∑​ck​≥∥z∥A​​(22)

结合(1)和(2)，我们得到：$F=\Vert \mathbf{z}{\Vert }_{\mathcal{A}}$。根据式(22)，只有当${\hat{p}}_k={\hat{c}}_k=|{\hat{s}}_k|$，且$\Vert \mathbf{z}{\Vert }_{\mathcal{A}}={\sum }_k{\hat{c}}_k={\sum }_k|{\hat{s}}_k|$时取等。因此，式(20)所示的原子分解可以直接得到原子范数。

优化问题(15)与优化问题(17)的关系：(17)中的SDP优化问题本质上是式(15)秩最小化问题的凸松弛。具体来讲，式(17)中目标函数的第二项$\frac{1}{2}{u}_1$是半正定Toeplitz矩阵$\mathbf{T}(\mathbf{u})$的核范数（所有奇异值之和）或者trace norm (up to a scaling factor)，which is a commonly used convex relaxation of the matrix rank。而第一项$\frac{1}{2}x$是为了防止平凡解(trivial solution)的出现设置的正则项。

关于问题求解过程的总结：
S1: the frequencies component are encoded in the Toeplitz matrix $\mathbf{T}(\mathbf{u})$
S2: Present the SDP optimization problem, and solve SDP. It follows that $\mathbf{u}$ is derived, i.e. $\mathbf{T}(\mathbf{u})$ is obtained.
S3: Once the resulting SDP problem is solved, the frequencies can be retrieved from the Vandermonde decomposition of $\mathbf{T}(\mathbf{u})$

另外，与原子${\ell }_0$范数类似，原子范数也可以被看作是基于协方差的方法，两者的不同之处在于原子${\ell }_0$范数直接利用了低秩性，而利用原子范数，问题更容易求解。

参考

[1] Yang, Z., Li, J., Stoica, P., & Xie, L. (2016). Sparse Methods for Direction-of-Arrival Estimation. ArXiv, abs/1609.09596.