Matrix pencil矩阵铅笔算法(原始论文记录与复现)（一）

《Estimating two-dimensional frequencies by matrix enhancement and matrix pencil》¹

---

问题概述

2D频率估计问题

从多径2D频率信号的$\sum$和式中分离出所有的频率分量的问题。

我们考虑一个无噪的二维频率数据域的2D-sample结构表达式——
$$x\left(m;n\right)=\sum _{i=1}^I{r}_i\exp \left(j{\varphi }_i+j2\pi f_{1i}m+j2\pi f_{2i}n\right)\left(1\right)$$
其中，$0⩽m⩽M,0⩽n⩽N.$
(1)表明2D频率的数据由I个多径的2D-frequencies参数：
{ ( f 1 i , f 2 i ) ; i = 1 , ⋯   , I } , { r i : i = 1 , ⋯   , I } , { ϕ i ; i = 1 , ⋯   , I } \left\{ \left( f_{1i},f_{2i} \right) ;i=1,\cdots ,I \right\} ,\left\{ r_i:i=1,\cdots ,I \right\} ,\left\{ \phi _i;i=1,\cdots ,I \right\} {(f1i​,f2i​);i=1,⋯,I},{ri​:i=1,⋯,I},{ϕi​;i=1,⋯,I}
分别为(非零的)模值、频率和相位。这些参数需要被估计。
在噪声系统中——
$$x^{\prime }\left(m;n\right)=x\left(m;n\right)+w\left(m;n\right)$$

其中 $w\left(m;n\right)$为2D的噪声序列（实际系统往往存在噪声）。

从式子(1)来看，$r_i,{\varphi }_i$可以很容易得到一旦 { ( f 1 i , f 2 i ) ; i = 1 , ⋯   , I } \left\{ \left( f_{1i},f_{2i} \right) ;i=1,\cdots ,I \right\} {(f1i​,f2i​);i=1,⋯,I}被精确估计之后因为根据指数乘法的性质$r_i\exp \left(j{\varphi }_i\right)$可以被线性地抽离出来。

问题的矩阵分解和简化

以exp为最小粒度重写式子(1)，可以得到——
$$x\left(m;n\right)=\sum _{i=1}^I{a}_i{y}_i^m{z}_i^n\left(2\right)$$
其中$\{\begin{array}{c}{y}_i=\exp \left(j2\pi f_{1i}\right)(3)\\ z_i=\exp \left(j2\pi f_{2i}\right)(4)\\ a_i=r_i\exp \left(j{\varphi }_i\right)(5)\end{array}。$对$y_i$和$z_i$的估计就对$f_{1i}$和$f_{2i}$的估计。
已知的抽样信号矩阵的形式如下——
$$\mathbf{X}=\left[\begin{array}{cccc}x\left(0;0\right)& x\left(0;1\right)& \cdots & x\left(0;N-1\right)\\ x\left(1;0\right)& x\left(1;1\right)& \cdots & x\left(1;N-1\right)\\ ⋮& ⋮& \cdots & ⋮\\ x\left(M-1;0\right)& x\left(M-1;1\right)& \cdots & x\left(M-1;N-1\right)\end{array}\right]\left(6\right)$$
由(3)在(7)的形式中可以得到——
$$\begin{gathered} \mathbf{X}=\mathbf{YAZ}\left(7\right) \\ \mathbf{Y}=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ y_1& y_2& \cdots & y_I\\ ⋮& ⋮& \cdots & ⋮\\ y_1^{M-1}& y_2^{M-1}& \cdots & y_I^{M-1}\end{array}\right]\left(8\right) \\ \mathbf{A}=\mathrm{diag}\left[a_1,a_2,\cdots ,a_I\right]\left(9\right) \\ \mathbf{Z}=\left[\begin{array}{cccc}1& z_1& \cdots & z_1^{N-1}\\ 1& z_2& \cdots & z_2^{N-1}\\ ⋮& ⋮& \cdots & ⋮\\ 1& z_I& \cdots & z_I^{N-1}\end{array}\right]\left(10\right) \end{gathered}$$
(这种二维形式的矩阵分解和乘法形式是常见的可以记住。)

具体的推导过程:
$$\begin{gathered} {\mathbf{X}}_{M\times N}=\underset{I\times I}{\underbrace{{\left[\begin{array}{cccc}1& 1& \cdots & 1\\ y_1& y_2& \cdots & y_I\\ ⋮& ⋮& \cdots & ⋮\\ y_1^{M-1}& y_2^{M-1}& \cdots & y_I^{M-1}\end{array}\right]}_{M\times I}\left[\begin{array}{cccc}{a}_1& 0& \cdots & 0\\ 0& a_2& \cdots & 0\\ ⋮& ⋮& \cdots & ⋮\\ 0& 0& \cdots & a_I\end{array}\right]}}{\left[\begin{array}{cccc}1& z_1& \cdots & z_1^{N-1}\\ 1& z_2& \cdots & z_2^{N-1}\\ ⋮& ⋮& \cdots & ⋮\\ 1& z_I& \cdots & z_I^{N-1}\end{array}\right]}_{I\times N} \\ =\left[\begin{array}{cccc}{a}_1{y}_1^0& a_2{y}_2^0& \cdots & a_I{y}_I^0\\ a_1{y}_1^1& a_2{y}_2^1& \cdots & a_I{y}_I^1\\ ⋮& ⋮& \cdots & ⋮\\ a_1{y}_1^{M-1}& a_2{y}_2^{M-1}& \cdots & a_I{y}_I^{M-1}\end{array}\right]\left[\begin{array}{cccc}{z}_1^0& z_1& \cdots & z_1^{N-1}\\ z_2^0& z_2& \cdots & z_2^{N-1}\\ ⋮& ⋮& \cdots & ⋮\\ z_I^0& z_I& \cdots & z_I^{N-1}\end{array}\right] \\ =\left[\begin{array}{cccc}{\sum }_{i=1}^I{a}_i{y}_i^0{z}_i^0& {\sum }_{i=1}^I{a}_i{y}_i^0{z}_i^1& \cdots & {\sum }_{i=1}^I{a}_i{y}_i^0{z}_i^{N-1}\\ {\sum }_{i=1}^I{a}_i{y}_i^1{z}_i^0& {\sum }_{i=1}^I{a}_i{y}_i^1{z}_i^1& \cdots & {\sum }_{i=1}^I{a}_i{y}_i^1{z}_i^{N-1}\\ ⋮& ⋮& \cdots & ⋮\\ {\sum }_{i=1}^I{a}_i{y}_i^{M-1}{z}_i^0& {\sum }_{i=1}^I{a}_i{y}_i^{M-1}{z}_i^1& \cdots & {\sum }_{i=1}^I{a}_i{y}_i^{M-1}{z}_i^{N-1}\end{array}\right] \\ =\left[\begin{array}{cccc}x\left(0;0\right)& x\left(0;1\right)& \cdots & x\left(0;N-1\right)\\ x\left(1;0\right)& x\left(1;1\right)& \cdots & x\left(1;N-1\right)\\ ⋮& ⋮& \cdots & ⋮\\ x\left(M-1;0\right)& x\left(M-1;1\right)& \cdots & x\left(M-1;N-1\right)\end{array}\right]\text{证毕}. \end{gathered}$$
·顺便复习一下考研中关于矩阵乘法$A_{m\times n}{B}_{n\times p}=C_{m\times p}$的行列分块和秩的结论。
对$A$做列分块，于是有
$$\begin{gathered} \left[{\mathbf{r}}_1,\cdots ,{\mathbf{r}}_n\right]\left[\begin{array}{ccc}{b}_{11}& \cdots & b_{1p}\\ ⋮& ⋮& ⋮\\ b_{n1}& \cdots & b_{np}\end{array}\right]=\left[{\mathbf{c}}_1,\cdots ,{\mathbf{c}}_p\right] \\ \{\begin{array}{c}{b}_{11}{\mathbf{r}}_1+b_{21}{\mathbf{r}}_2+\cdots +b_{n1}{\mathbf{r}}_n={\mathbf{c}}_1\\ ⋮\\ b_{1i}{\mathbf{r}}_1+b_{2i}{\mathbf{r}}_2+\cdots +b_{ni}{\mathbf{r}}_n={\mathbf{c}}_i\\ ⋮\\ b_{1p}{\mathbf{r}}_1+b_{2p}{\mathbf{r}}_2+\cdots +b_{np}{\mathbf{r}}_n={\mathbf{c}}_p\end{array} \end{gathered}$$
$C=AB$的列向量可以由$A$的列向量线性表出。同时表明如果$A$列满秩，$r(AB)=r(B)$。即左乘列满秩矩阵秩不变。
对$B$做行分块，于是有
$$\begin{gathered} \left[\begin{array}{ccc}{a}_{11}& \cdots & a_{1n}\\ ⋮& ⋮& ⋮\\ a_{m1}& \cdots & a_{mn}\end{array}\right]\left[\begin{array}{c}{\mathbf{b}}_1\\ ⋮\\ {\mathbf{b}}_n\end{array}\right]=\left[\begin{array}{c}{c}_1\\ ⋮\\ c_m\end{array}\right] \\ \{\begin{array}{c}{a}_{11}{\mathbf{b}}_1+a_{12}{\mathbf{b}}_2+\cdots +a_{1n}{\mathbf{b}}_n={\mathbf{c}}_1\\ ⋮\\ a_{i1}{\mathbf{b}}_1+a_{i2}{\mathbf{b}}_2+\cdots +a_{in}{\mathbf{b}}_n={\mathbf{c}}_i\\ ⋮\\ a_{m1}{\mathbf{b}}_1+a_{m2}{\mathbf{b}}_2+\cdots +a_{mn}{\mathbf{b}}_n={\mathbf{c}}_m\end{array} \end{gathered}$$
$C=AB$的行向量可以由$B$的列向量线性表出。同时表明如果$B$行满秩，$r(AB)=r(A)$。即右乘行满秩矩阵秩不变。
同时表明，对于式子(3)矩阵$X_{m,n}={\sum }_{i=1}^I{a}_i{y}_i^m{z}_i^n$来说，$z_i^n$ 是来自于(10)中$Z$的第$n$列的贡献，$y_i^m$ 是来自于(10)中$Y$的第$m$行的贡献，对角矩阵$A$则是将相应的$y_i^m$和$z_i^n$ $concat$起来，于是上面式子的由来被完整地梳理了一遍。

式子(7)$⟶$ $r\left(YAZ\right)\ll r\left(A\right)=I$ ,等号当且仅当$Y$列满秩、$Z$行满秩且$M⩾I,N⩾I$. 观察到$Y$和$Z$矩阵的$VanderMonde$形式，上述条件可以被恒等转化为 { y i ; i = 1 , ⋯   , I } \left\{ y_i;i=1,\cdots ,I \right\} {yi​;i=1,⋯,I}和 { z i ; i = 1 , ⋯   , I } \left\{ z_i;i=1,\cdots ,I \right\} {zi​;i=1,⋯,I}都不包含重复值否则对于$r(X)=I$这个条件的推理的不够充分的(是一种病态条件)。ps:原文中作者认为只要2D频率对 { { y i , z i } ; i = 1 , ⋯   , I } \left\{ \left\{ y_i,z_i \right\} ;i=1,\cdots ,I \right\} {{yi​,zi​};i=1,⋯,I}的每组的独立的就ok，也就是允许$y$和$z$的单独频率可以重复但需要错开，我认为这种方法不能够保证Y和Z的满秩条件，在仿真的实验中也验证了$y$和$z$都需要分别包含不重复值这一结论，评论区有对这个问题思考的可以一起交流！
(update:这里后来重新理解了一下就是式子(7)要使$X$秩为$I$确实要y和z都不包含相同值，但是从后面增强矩阵式子(25)可以看出只要yz对没有重复值就可以很好地恢复，所以在这个意义上 这种增强强化了MEMP算法的使用范畴比较2D对重复的可能性比1D频率重复的可能性小得多。)

一、算法推导与展开

matrix enhancement

增强型矩阵${\mathbf{X}}_e$

通过矩阵的分隔和合并来构造增强后的Hankle-block矩阵(每一条副对角线上的元素都相等)
$${\mathbf{X}}_e={\left[\begin{array}{cccc}{\mathbf{X}}_0& {\mathbf{X}}_1& \cdots & {\mathbf{X}}_{M-K}\\ {\mathbf{X}}_1& {\mathbf{X}}_2& \cdots & {\mathbf{X}}_{M-K+1}\\ ⋮& ⋮& \cdots & ⋮\\ {\mathbf{X}}_{K-1}& {\mathbf{X}}_K& \cdots & {\mathbf{X}}_{M-1}\end{array}\right]}_{KL\times \left(M-K+1\right)\left(N-L+1\right)}\left(11\right)$$
其中——
$${\mathbf{X}}_m={\left[\begin{array}{cccc}x\left(m;0\right)& x\left(m;1\right)& \cdots & x\left(m;N-L\right)\\ x\left(m;1\right)& x\left(m;2\right)& \cdots & x\left(m;N-L+1\right)\\ ⋮& ⋮& \cdots & ⋮\\ x\left(m;L-1\right)& x\left(m;L\right)& \cdots & x\left(m;N-1\right)\end{array}\right]}_{L\times \left(N-L+1\right)}\left(12\right)$$
${\mathbf{X}}_m$是给定$m$之后对$x\left(m;0\right)⟶x\left(m;N-1\right)$以$L$为长度做滑动的矩阵。
${\mathbf{X}}_e$是$m$取值在$0⟶M-1$之间的以$K$为长度做滑动的${\mathbf{X}}_m$。
式子(2)带入(12)，
$$\begin{gathered} {\mathbf{X}}_m={\mathbf{Z}}_L{\mathbf{AY}}_d^m{\mathbf{Z}}_R\left(13\right) \\ {\mathbf{Z}}_L={\left[\begin{array}{cccc}{z}_1^0& z_2^0& \cdots & z_I^0\\ z_1& z_2& \cdots & z_I\\ ⋮& ⋮& \cdots & ⋮\\ z_1^{L-1}& z_2^{L-1}& \cdots & z_I^{L-1}\end{array}\right]}_{L\times I}\left(14\right) \\ {\mathbf{Y}}_d=\mathrm{diag}\left[y_1,y_2,\cdots ,y_I\right]\left(15\right) \\ {\mathbf{Z}}_R={\left[\begin{array}{cccc}{z}_1^0& z_1& \cdots & z_1^{N-L}\\ z_2^0& z_2& \cdots & z_2^{N-L}\\ ⋮& ⋮& \cdots & ⋮\\ z_I^0& z_I& \cdots & z_I^{N-L}\end{array}\right]}_{I\times N-L+1}\left(16\right) \end{gathered}$$
式子(13)的由来和(7)类似，对于式子(12)的矩阵${{\mathbf{X}}_m}_{\left(pq\right)}=x\left(m;p+q\right)={\sum }_{i=1}^I{a}_i{y}_i^m{z}_i^{p+q}$, 其中$p,q$从0开始索引，$a_i{y}_i^m$由两个对角矩阵的乘积${\mathbf{AY}}_d^m$贡献，$z_i^{p+q}$由第$p$行的${\mathbf{Z}}_L$和第$q$列的${\mathbf{Z}}_R$的内积贡献，以上思路同式子(7)的由来。
将(13)带入(11),
$$\begin{gathered} {\mathbf{X}}_e={\mathbf{E}}_L{\mathbf{AE}}_R\left(17\right) \\ {\mathbf{E}}_L={\left[\begin{array}{l}{\mathbf{Z}}_L\\ {\mathbf{Z}}_L{\mathbf{Y}}_d\\ \cdots \\ {\mathbf{Z}}_L{\mathbf{Y}}_d^{K-1}\end{array}\right]}_{KL\times I}\left(18\right) \\ {\mathbf{E}}_R={\left[{\mathbf{Z}}_R,{\mathbf{Y}}_d{\mathbf{Z}}_R,\cdots ,{\mathbf{Y}}_d^{M-K}{\mathbf{Z}}_R\right]}_{I\times \left(N-L+1\right)\left(M-K+1\right)}\left(19\right) \end{gathered}$$
式子(17)的由来还是类似上面，
$${{\mathbf{X}}_e}_{\left(pq\right)}={\mathbf{X}}_{p+q-2}={\mathbf{Z}}_L{\mathbf{AY}}_d^{p+q-2}{\mathbf{Z}}_R={\mathbf{Z}}_L{\mathbf{Y}}_d^{p-1}{\mathbf{AY}}_d^{q-1}{\mathbf{Z}}_R={\mathbf{E}}_L\left(p,:\right)\cdot \mathbf{A}\cdot {\mathbf{E}}_R\left(:,q\right)\left(20\right)$$
式子(17)表明增强后的矩阵${\mathbf{X}}_e$满足
$$r\left({\mathbf{X}}_e\right)=I\iff r\left({\mathbf{E}}_L\right)=r\left({\mathbf{E}}_R\right)=I$$

$K$和$L$满足的条件

$$Theorm1:r\left({\mathbf{E}}_L\right)=I\iff K⩾IandL\gg I.$$

原文中引入了一个shuffle重拍列矩阵
$$\mathbf{P}={\left[\begin{array}{c}\begin{array}{c}{\mathbf{p}}^T\left(1\right)\\ {\mathbf{p}}^T\left(1+L\right)\\ \cdots \\ {\mathbf{p}}^T\left(1+\left(K-1\right)L\right)\end{array}\\ \begin{array}{c}{\mathbf{p}}^T\left(2\right)\\ {\mathbf{p}}^T\left(2+L\right)\\ \cdots \\ {\mathbf{p}}^T\left(2+\left(K-1\right)L\right)\\ \cdots \\ \cdots \\ \begin{array}{c}{\mathbf{p}}^T\left(L\right)\\ {\mathbf{p}}^T\left(L+L\right)\\ \cdots \\ {\mathbf{p}}^T\left(L+\left(K-1\right)L\right)\end{array}\end{array}\end{array}\right]}_{KL\times KL}\left(21\right)$$
其中$\mathbf{p}\left(i\right)$表示一个长度为$KL$的列向量，其中位置是$i$为1，其余为0,易得
$$\begin{gathered} {\mathbf{E}}_{LP}={\mathbf{PE}}_L\left(22\right) \\ {\mathbf{E}}_{LP}={\left[\begin{array}{c}{\mathbf{Y}}_L\\ {\mathbf{Y}}_L{\mathbf{Z}}_d\\ \cdots \\ {\mathbf{Y}}_L{\mathbf{Z}}_d^{L-1}\end{array}\right]}_{KL\times I}\left(23\right) \\ {\mathbf{Y}}_L={\left[\begin{array}{cccc}1& 1& \cdots & 1\\ y_1& y_2& \cdots & y_I\\ ⋮& ⋮& \cdots & ⋮\\ y_1^{L-1}& y_2^{L-1}& \cdots & y_I^{L-1}\end{array}\right]}_{L\times I}\left(24\right) \\ {\mathbf{Z}}_d=\mathrm{diag}\left(z_1,z_2,\cdots ,z_I\right)\left(25\right) \end{gathered}$$
$\mathbf{P}$的作用就是交换${\mathbf{E}}_L$中y和z的位置，即$y_i$在${\mathbf{E}}_L$的位置相当于$z_i$在${\mathbf{E}}_{LP}$. 这种重拍变换的逻辑在于在${\mathbf{E}}_L$的矩阵结构中——
$$for{\mathbf{E}}_L={\left[\begin{array}{l}{\mathbf{Z}}_L\\ {\mathbf{Z}}_L{\mathbf{Y}}_d\\ \cdots \\ {\mathbf{Z}}_L{\mathbf{Y}}_d^{L-1}\end{array}\right]}_{KL\times I}=\left[\begin{array}{l}\left[\begin{array}{cccc}{y}_1^0& y_2^0& \cdots & y_I^0\\ z_1{y}_1^0& z_2{y}_2^0& \cdots & z_I{y}_I^0\\ ⋮& ⋮& ⋮& ⋮\\ z_1^{L-1}{y}_1^0& z_2^{L-1}{y}_2^0& \cdots & z_I^{L-1}{y}_I^0\end{array}\right]\\ \left[\begin{array}{cccc}{y}_1^1& y_2^1& \cdots & y_I^1\\ z_1{y}_1^1& z_2{y}_2^1& \cdots & z_I{y}_I^1\\ ⋮& ⋮& ⋮& ⋮\\ z_1^{L-1}{y}_1^1& z_2^{L-1}{y}_2^1& \cdots & z_I^{L-1}{y}_I^1\end{array}\right]\\ \cdots \\ \left[\begin{array}{cccc}{y}_1^{K-1}& y_2^{K-1}& \cdots & y_I^{K-1}\\ z_1{y}_1^{K-1}& z_2{y}_2^{K-1}& \cdots & z_I{y}_I^{K-1}\\ ⋮& ⋮& ⋮& ⋮\\ z_1^{L-1}{y}_1^{K-1}& z_2^{L-1}{y}_2^{K-1}& \cdots & z_I^{L-1}{y}_I^{K-1}\end{array}\right]\end{array}\right]\left(25\right)$$
相当于取每个矩阵的第一行拼接在一起，然后第二行，以此类推，相当于对行进行了重排。这样就把y和z的位置交换，同时也把K和L交换。
这样${\mathbf{Z}}_L$和${\mathbf{Y}}_L$都是${\mathbf{E}}_L$的子矩阵，于是有
$$r\left({\mathbf{E}}_L\right)⩾r\left[\begin{array}{c}{\mathbf{Z}}_L\\ {\mathbf{Y}}_L\end{array}\right]=r\left[\begin{array}{c}{\left[\begin{array}{cccc}{z}_1^0& z_2^0& \cdots & z_I^0\\ z_1& z_2& \cdots & z_I\\ ⋮& ⋮& \cdots & ⋮\\ z_1^{L-1}& z_2^{L-1}& \cdots & z_I^{L-1}\end{array}\right]}_{L\times I}\\ {\left[\begin{array}{cccc}{y}_1^0& y_2^0& \cdots & y_I^0\\ y_1& y_2& \cdots & y_I\\ ⋮& ⋮& \cdots & ⋮\\ y_1^{L-1}& y_2^{L-1}& \cdots & y_I^{L-1}\end{array}\right]}_{L\times I}\end{array}\right]\left(26\right)$$
如果 { ( y i , z i ) ; i = 1 , 2 , ⋯   , I } i s    d i s t i n c t \left\{ \left( y_i,z_i \right) ;i=1,2,\cdots ,I \right\} is\,\,distinct {(yi​,zi​);i=1,2,⋯,I}isdistinct. (26)矩阵列满秩
$\therefore$当$K⩾IandL\gg I⟶r\left({\mathbf{E}}_L\right)=I$.充分性得证。
关于必要性，由于${\mathbf{E}}_L$是$KL$行$I$列矩阵，所以$r\left({\mathbf{E}}_L\right)=I$必须有$KL⩾I$，这个条件不够强，可能会失效，而充分条件是最严格的，即
$$\begin{gathered} K⩾IandL\gg I\left(27\right) \\ KL\gg I\left(28\right) \end{gathered}$$
根据对称性可以得到${\mathbf{E}}_R$侧的结论，
$$\begin{gathered} (M-K+1)⩾Iand(N-L+1)\gg I\left(29\right) \\ (M-K+1)(N-L+1)\gg I\left(30\right) \end{gathered}$$
事实上, ${\mathbf{E}}_L$和${\mathbf{E}}_R$是行拼接和列拼接分别对应列满秩和行满秩的特征，${\mathbf{E}}_R$侧的$\mathbf{P}$矩阵是做列shuffle重排！

${\mathbf{X}}_e$的特征结构

我们考虑从${\mathbf{X}}_e$中提取2D的频率信息——
对${\mathbf{X}}_e$做SVD分解，得到

$$\begin{gathered} {\mathbf{X}}_e=\sum _{i=1}^{\min }{\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^H={\mathbf{U}}_s{\Sigma }_s{\mathbf{V}}_s^H+{\mathbf{U}}_n{\Sigma }_n{\mathbf{V}}_n^H\left(31\right) \\ {\mathbf{U}}_s=\left[{\mathbf{u}}_1,{\mathbf{u}}_2,\cdots ,{\mathbf{u}}_I\right]\left(31.a\right) \\ {\Sigma }_s=\mathrm{diag}\left[{\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_I\right]\left(31.b\right) \\ {\mathbf{V}}_s=\left[{\mathbf{v}}_1,{\mathbf{v}}_2,\cdots ,{\mathbf{v}}_I\right]\left(31.c\right) \\ {\mathbf{U}}_n=\left[{\mathbf{u}}_{I+1},{\mathbf{u}}_{I+2},\cdots ,{\mathbf{u}}_{\min }\right]\left(31.d\right) \\ {\Sigma }_n=\mathrm{diag}\left[{\sigma }_{I+1},{\sigma }_{I+2},\cdots ,{\sigma }_{\min }\right]\left(31.e\right) \\ {\mathbf{V}}_n=\left[{\mathbf{v}}_{I+1},{\mathbf{v}}_{I+2},\cdots ,{\mathbf{v}}_{\min }\right]\left(31.f\right) \end{gathered}$$
这里${\sigma }_1⩾{\sigma }_2⩾\cdots ⩾{\sigma }_{\min }$.对于无噪声的情况，只有$I$个非零奇异值，而有噪声情况下，噪声带来的奇异值的扰动会和2D频率的奇异值相差较大从而可以分离出多径，从而预测$I$。
由于$r\left(X_e\right)=r\left(E_L\right)=I$,所以
$$\begin{gathered} range\left(X_e\right)=range\left(E_L\right)=range\left(U_s\right) \\ range\left(X_e^H\right)=range\left(E_R^H\right)=range\left(V_s\right)\left(32\right) \end{gathered}$$
这是因为式子(17)中${\mathbf{E}}_L$列满秩，${\mathbf{E}}_R$行满秩，所以SVD的$U$矩阵的列空间和原矩阵的列空间一致，$V$矩阵的行空间和原矩阵的行空间一致。
再次考虑(18)(25)中${\mathbf{E}}_L$的结构, ${\mathbf{E}}_L$的第$i$列可以被表示为
$${\mathbf{e}}_{Li}={\mathbf{y}}_{Li}\otimes {\mathbf{z}}_{Li}\left(33\right)$$
其中$\otimes$代表Kronecker product。
这代表${\mathbf{E}}_L$的列向量可以和$y_i$和$z_i$的克罗内克积。式子(33)表明了${\mathbf{E}}_L$、${\mathbf{y}}_L$和${\mathbf{z}}_L$之间的关系，
$${\mathbf{e}}_L={\mathbf{y}}_L\otimes {\mathbf{z}}_L\left(34\right)$$
其中，
$$\begin{gathered} {\mathbf{y}}_L={\left[1,y,\cdots ,y^{K-1}\right]}^T\left(35\right) \\ {\mathbf{z}}_L={\left[1,z,\cdots ,z^{L-1}\right]}^T\left(36\right) \\ y=\exp \left(j2\pi f_1\right)\left(37\right) \\ z=\exp \left(j2\pi f_2\right)\left(38\right) \end{gathered}$$
可以发现
e L ∈ { e L 1 , e L 2 , ⋯   , e L I } ⇔ ( f 1 , f 2 ) = ( f 1 i , f 2 i ) \boldsymbol{e}_L\in \left\{ \boldsymbol{e}_{L1},\boldsymbol{e}_{L2},\cdots ,\boldsymbol{e}_{LI} \right\} \Leftrightarrow \left( f_1,f_2 \right) =\left( f_{1i},f_{2i} \right) eL​∈{eL1​,eL2​,⋯,eLI​}⇔(f1​,f2​)=(f1i​,f2i​)
因为(25)${\mathbf{E}}_L$的结构天然地将对应的$y_i$和$z_i$联系在一起。
另外一个很重要的信息就是y和z正确配对后有${\mathbf{e}}_L\perp {\mathbf{U}}_n$,所以只需要最大化下面的式子就可以找到对应的2D频率——
$$Maximize⟶\frac{1}{{\sum }_{i=I+1}^{\min }||{\mathbf{u}}_i^H{\mathbf{e}}_L\left(f_1,f_2\right)||}\left(39\right)$$

matrix pencil

Matrix pencil在Gene H.Golub教授的《Matrix Computations》²中的描述如下：

If $A,B\in {\mathbb{C}}^{n\times n}$，then the set of all matrices of the form $A-\lambda B$ with $\lambda \in \mathbb{C}$ is a $pencil$.
The $generalizedeigenvalues$ of $A-\lambda B$ are elements of the set $\lambda \left(A,B\right)$ defined by
λ ( A , B ) = { z ∈ C : det ⁡ ( A − z B ) = 0 } \lambda \left( A,B \right) =\left\{ z\in \mathbb{C} :\det \left( A-zB \right) =0 \right\} λ(A,B)={z∈C:det(A−zB)=0}
If $\lambda \in \lambda \left(A,B\right)and0\ne x\in {\mathbb{C}}^{n}satisfies$.
$$Ax=\lambda Bx$$
Then $x$ is an $eigenvector$ of $A-\lambda B$.

矩阵铅笔方法可以看成将两个矩阵以某种方式重构使得需要估计的参数等于矩阵铅笔的秩的减少数(如广义特征值)。

提取$y_i$

假设充分条件(29)可以被满足，则$range\left({\mathbf{U}}_s\right)=range\left({\mathbf{E}}_L\right)$，于是
有且仅有唯一一个$I\times I$的非奇异矩阵$\mathbf{T}$,使得${\mathbf{U}}_s={\mathbf{E}}_L\mathbf{T}.$ $\mathbf{T}$的唯一性可以由${\mathbf{U}}_s$与${\mathbf{E}}_L$列满秩(均包含独立的$I$列)得到。
$$\{\begin{array}{c}{\mathbf{U}}_1={\mathbf{U}}_s\text{但是不包含}{U}_s\text{的最后}L\text{行}\\ {\mathbf{U}}_2={\mathbf{U}}_s\text{但是不包含}{U}_s\text{的开头}L\text{行}\end{array}\{\begin{array}{c}{\mathbf{U}}_1={\mathbf{E}}_1\mathbf{T}\\ {\mathbf{U}}_2={\mathbf{E}}_1{\mathbf{Y}}_d\mathbf{T}\end{array}$$
${\mathbf{U}}_2$和${\mathbf{U}}_1$构成的矩阵铅笔包含了
$${\mathbf{U}}_2-\lambda {\mathbf{U}}_1={\mathbf{E}}_1\left({\mathbf{Y}}_d-\lambda \mathbf{I}\right)\mathbf{T}$$
所以满足${\mathbf{U}}_2-\lambda {\mathbf{U}}_1=0$的广义特征值就是
$${\mathbf{Y}}_d=\mathrm{diag}\left[y_1,y_2,\cdots ,y_I\right]。$$
也就是说$y_i$就是矩阵铅笔${\mathbf{U}}_2-\lambda {\mathbf{U}}_1$减少的秩(因为当${\lambda }_i=y_i$矩阵铅笔的秩就会减一)。上述推理要求矩阵${\mathbf{E}}_1$和$\mathbf{T}$均列满秩为$I$.

提取$z_i$

矩阵$\mathbf{P}$将${\mathbf{E}}_L$重排成${\mathbf{E}}_{LP}$，于是将 ${\mathbf{U}}_s$重排成${\mathbf{U}}_{sP}$
$$\begin{gathered} {\mathbf{U}}_{sP}={\mathbf{PU}}_s \\ \{\begin{array}{c}{\mathbf{U}}_{1P}={\mathbf{U}}_{sP}\text{但是不包含}{\mathbf{U}}_{sP}\text{的最后}K\text{行}\\ {\mathbf{U}}_{2P}={\mathbf{U}}_{sP}\text{但是不包含}{\mathbf{U}}_{sP}\text{的开头}K\text{行}\end{array} \\ {\mathbf{U}}_{2P}-\lambda {\mathbf{U}}_{1P}={\mathbf{E}}_{1P}\left({\mathbf{Z}}_d-\lambda \mathbf{I}\right)\mathbf{T} \end{gathered}$$
所以满足${\mathbf{U}}_{2P}-\lambda {\mathbf{U}}_{1P}=0$的广义特征值就是

$${\mathbf{Z}}_d=\mathrm{diag}\left[z_1,z_2,\cdots ,z_I\right].$$

$K$和$L$满足的条件

$$\{\begin{array}{c}r\left({\mathbf{E}}_1\right)=I⟶\{\begin{array}{c}K-1⩾I\\ L⩾I\end{array}\\ r\left({\mathbf{E}}_{1P}\right)=I⟶\{\begin{array}{c}K⩾I\\ L-1⩾I\end{array}\end{array}\left(40\right)$$
结合(40)和(29)得到MEMP的总体充分条件：
$$\{\begin{array}{c}M-I+1>K⩾I+1\\ N-I+1⩾L-I+1\end{array}\left(41\right)$$

必要条件是
$$\{\begin{array}{c}\left(K-1\right)L⩾I\\ K\left(L-1\right)⩾I\\ \left(M-K+1\right)\left(N-L+1\right)⩾I\end{array}\left(42\right)$$

2D 匹配pairing

现在可以通过matrix pencil构造出2D频率就是${\mathbf{U}}_2-\lambda {\mathbf{U}}_1$和
${\mathbf{U}}_{2P}-\lambda {\mathbf{U}}_{1P}$的广义特征值，但不一定是对应匹配的。
(39)表明，对于$I$组2D频率pair，我们应该要

$$Minimize⟶J_n\left(i,j\right)=\sum _{i=I+1}^{\min }||{\mathbf{u}}_i^H{\mathbf{e}}_L\left(y_i,z_j\right)|{|}^2\left(43\right)$$
由SVD特性，
$$Maximize⟶J_n\left(i,j\right)=\sum _{i=1}^I||{\mathbf{u}}_i^H{\mathbf{e}}_L\left(y_i,z_j\right)|{|}^2\left(44\right)$$
所以遍历$y_i$和$z_i$使得(44)满足即可。

---

总结

---

1. Y. Hua, “Estimating two-dimensional frequencies by matrix enhancement and matrix pencil,” in IEEE Transactions on Signal Processing, vol. 40, no. 9, pp. 2267-2280, Sept. 1992, doi: 10.1109/78.157226. ↩︎

2. Matrix Computations - 4th Edition. Gene H.Golub ↩︎