PAST 算法详细推导

本文源于《矩阵分析与应用》的这样一句话：

如何运用矩阵求逆引理得到 PAST 算法，书中省略了推导过程。今尝试推导如下。

考虑特征子空间的跟踪与更新问题，「投影逼近」是解决此问题的一个方法。将特征子空间的确定当作一个无约束最优化问题来求解，称为「投影逼近子空间跟踪」（Projection Approximation Subspace Tracking，PAST）。

特征子空间的跟踪与更新问题转化为以 J(W) 为目标函数的极小化问题

\[ J(W) = E \left\{ \| x - WW^Hx \|^2 \right\} \]

其中，\(x\) 是 \(n \times 1\) 观测数据向量，\(W\) 是 \(n \times r\) 矩阵。求矩阵 \(W\) 使得目标函数达到全局极小点 \(\text{min}J(W)\) 。

已被证明，当 \(W\) 的列向量张成的子空间等于信号子空间时，目标函数 \(J(W)\) 达到全局极小点。所谓「子空间跟踪」其跟踪的就是能够使 \(W\) 成为信号子空间的 \(W\)。

什么是信号子空间

何种情况下 \(W\) 会成为信号子空间？这与 PAST 算法本身无关，而属于子空间分析方法中的内容。在子空间分析中，观测数据向量 \(x\) 的自相关矩阵表示为 \(R = E\{xx^H\}\)，对相关矩阵进行特征值分解 \(R = U \Sigma U^H\)。矩阵 \(U\) 的列向量分为两个部分 \(U = [U_S, U_G] = [u_1, u_2, ..., u_r, u_{r+1}, ..., u_n]\)。自相关矩阵的 \(r\) 个主特征向量 \(U_S = [u_1, u_2, ..., u_r]\) 张成的子空间即为观测数据空间的信号子空间，剩余的 \(n-r\) 个次特征向量 \(U_G = [u_{r+1}, ..., u_n]\) 张成的子空间即为观测数据空间的噪声子空间。

信号子空间与投影矩阵

信号子空间与噪声子空间是正交的，即有 \(\text{Span}\{u_1, u_2, ..., u_r\} \perp \text{Span}\{u_{r+1}, ..., u_n\}\)。由于 \(U\) 是酉矩阵，因此有

\[ UU^H = U_S U_S^H + U_G U_G^H = I \]

定义信号子空间上的投影矩阵为

\[ P_S = U_S (U_S^H U_S)^{-1} U_S^H \] 只有当满足正交性约束条件 \(U_S^H U_S = I\) 时，信号子空间上的投影矩阵可简化为 \(P_S = U_S U_S^H\)。于是，\(P_Sx = U_S U_S^Hx\) 为向量 \(x\) 在信号子空间上的投影。

同理，定义噪声子空间上的投影矩阵为

\[ P_G = U_G (U_G^H U_G)^{-1} U_G^H \] 只有当满足正交性约束条件 \(U_G^H U_G = I\) 时，噪声子空间上的投影矩阵可简化为 \(P_G = U_G U_G^H\)。于是，\(P_Gx = U_G U_G^Hx\) 为向量 \(x\) 在噪声子空间上的投影。

由于信号子空间与噪声子空间具有正交关系，因此 \(P_S\) 和 \(P_G\) 互为正交投影矩阵，即 \(P_S + P_G = I\)。因此 \(P_G = (I - P_S)x\) 还表示向量 \(x\) 在信号子空间上的正交投影。

另有如下子空间分析的引理，其中 \(P_{H_1}\) 是到子空间 \(H_1\) 的投影矩阵，\(P_{H_2}\) 是到子空间 \(H_2\) 的投影矩阵。

此引理表明，若一个矩阵等于另一个矩阵右乘一个正交矩阵，则这两个矩阵等价或张成相同的列空间。这就引出了第二个约束条件，即齐次性约束条件。

极小化问题解的形式

对于目标函数 \(J(W)\)，若矩阵 \(W\) 满足 1）正交性约束条件 \(W^H W = I\)，2）齐次性约束条件 \(J(W) = J(WQ)\)，其中 \(Q\) 为 \(r \times r\) 任意正交矩阵，则极小化问题 \(\text{min} J(W)\) 的解不是一个 \(W\) 矩阵，而是由 \(WQ\) 组成的矩阵集合。

从子空间分析到 PAST 算法

PAST 算法给出了求解极小化问题 \(\text{min} J(W)\) 的方法。已被证明，当且仅当 \(W = U_r Q\)，其中 \(Q\) 是 \(r \times r\) 任意酉矩阵，而 \(U_r\) 是 \(n \times r\) 矩阵，其 \(r\) 个列向量是自相关矩阵 \(R = E\{xx^H\}\) 的 \(r\) 个主特征向量，则 \(W\) 的列向量张成的子空间为信号子空间，此时目标函数 \(J(W)\) 达到全局极小值。

更重要的是，这一结论的得出并不事先假定有任何约束条件，而是在目标函数 \(J(W)\) 达到极小值时自动满足 1）正交性约束条件 \(W^H W = I\)，2）齐次性约束条件 \(J(W) = J(WQ)\)。这就是为什么说，原本具有正交性约束和齐次性的极小化问题可以等价为一个无约束的最优化问题。如此一来，求极小化问题 \(\text{min} J(W)\) 的解即转化为求自相关矩阵的 \(r\) 个主特征向量。

然而在实际应用中，自相关矩阵是随时间变化的，其特征值和特征向量也是随时间变化的。跟踪时变信号子空间 \(W\) 的一个办法是采用递推最小二乘法。

PAST 算法推导

定义指数加权的目标函数 \[ \begin{equation} \tag{1} \begin{aligned} J_1(W(t)) &= \sum_{i = 1}^{t} \beta^{t - i} \|x(i) - W(t)W^H(t)x(i)\|^2 \\ &= \sum_{i = 1}^{t} \beta^{t - i} \|x(i) - W(t)y(i)\|^2 \end{aligned} \end{equation} \]

其中，\(y(i) = W(t) x(i)\)。

由自适应滤波理论知，极小化问题 \(\text{min} J_1(W)\) 的最优解为维纳滤波器。

\[ W(t) = C_{xy}(t) C_{yy}^{-1}(t) \tag{2} \]

其中，互相关矩阵 \(C_{xy}(t)\) 和自相关矩阵 \(C_{yy}(t)\) 可以递推得到

\[ C_{xy}(t) = \sum_{i = 1}^{t} \beta^{t - i} x(i) y^H(t) = \beta C_{xy}(t - 1) + x(t) y^H(t) \tag{3} \]

\[ C_{yy}(t) = \sum_{i = 1}^{t} \beta^{t - i} y(i) y^H(t) = \beta C_{yy}(t - 1) + y(t) y^H(t) \tag{4} \]

对式（4）应用矩阵求逆引理，令

\[ \begin{equation} \tag{5} \begin{aligned} A &= C_{yy}(t) \\ B^{-1} &= \beta C_{yy}(t - 1) \\ C &= y(t) \\ D &= 1 \end{aligned} \end{equation} \]

则相关矩阵的逆矩阵为

\[ C_{yy}^{-1}(t) = \beta^{-1} C_{yy}^{-1}(t - 1) - \frac{\beta^{-1} C_{yy}^{-1}(t - 1) y(t) · \beta^{-1} y^H(t) C_{yy}^{-1}(t - 1)}{1 + \beta^{-1} y^H(t) C_{yy}^{-1}(t - 1) y(t)} \tag{6} \]

令 \(P(t) = C_{yy}^{-1}(t)\)，则上式可重写为

\[ P(t) = \beta^{-1} P(t - 1) - \frac{\beta^{-1} P(t - 1) y(t)}{1 + \beta^{-1} y^H(t) P(t - 1) y(t)} · \beta^{-1} y^H(t) P(t - 1) \tag{7} \]

令

\[ g(t) = \frac{\beta^{-1} P(t - 1) y(t)}{1 + \beta^{-1} y^H(t) P(t - 1) y(t)} \tag{8} \]

带入式（7）可得

\[ P(t) = \beta^{-1} P(t - 1) - \beta^{-1} g(t) y^H(t) P(t-1) \tag{9} \]

另外整理式（8）可得

\[ \begin{equation} \tag{10} \begin{aligned} g(t) &= \beta^{-1} P(t - 1) y(t) - \beta^{-1} g(t) y^H(t) P(t - 1) y(t) \\ &= [\beta^{-1} P(t - 1) - \beta^{-1} g(t) y^H(t) P(t - 1)] y(t) \\ &= P(t) y(t) \end{aligned} \end{equation} \]

式（8~10）导出了增益向量 \(g(t) = P(t) y(t)\)，后面将会再次用到它。现继续整理式（7）得

\[ P(t) = \beta^{-1} P(t - 1) - \frac{P(t - 1) y(t)}{\beta + y^H(t) P(t - 1) y(t)} · \beta^{-1} y^H(t) P(t - 1) \tag{11} \]

令 \(h(t) = P(t - 1) y(t)\)，由于自相关矩阵是 Hermitian 矩阵，亦有 \(h^H(t) = y^H(t) P(t - 1)\)，则上式可重写为

\[ P(t) = \beta^{-1} P(t - 1) - \beta^{-1} \frac{h(t)}{\beta + y^H(t) h(t)} h^H(t) \tag{12} \]

令 \(\alpha(t) = \beta + y^H(t) h(t)\)，则式（12）可重写为

\[ P(t) = \beta^{-1} [P(t - 1) - h(t) h^H(t) / \alpha(t)] \tag{13} \]

式（11 ~ 13）导出了自相关矩阵的逆矩阵 \(P(t)\) 的计算式。

下面使用 \(P(t) = C_{yy}^{-1}(t)\) 和式（3）将式（2）重写为

\[ \begin{equation} \tag{14} \begin{aligned} W(t) &= C_{xy}(t) C_{yy}^{-1}(t) \\ &= C_{xy}(t) P(t) \\ &= [\beta C_{xy}(t - 1) + x(t) y^H(t)] P(t) \\ &= \beta C_{xy}(t -1) P(t) + x(t) y^H(t) P(t) \end{aligned} \end{equation} \] 将式（13）带入上式右边第一项得

\[ \begin{equation} \tag{15} \begin{aligned} W(t) &= \beta C_{xy}(t -1)\{\beta^{-1} [P(t - 1) - h(t) h^H(t) / \alpha(t)]\} + x(t) y^H(t) P(t)\\ &= C_{xy}(t - 1) P(t - 1) - C_{xy}(t - 1) h(t) h^H(t) / \alpha(t) + x(t) y^H(t) P(t) \\ &= C_{xy}(t - 1) P(t - 1) - C_{xy}(t - 1) P(t - 1) y(t) h^H(t) / \alpha(t) + x(t) y^H(t) P(t) \\ &= W(t -1) - W(t -1) y(t) h^H(t) / \alpha(t) + x(t) y^H(t) P(t) \end{aligned} \end{equation} \]

前面由式（10）导出增益向量 \(g(t) = P(t) y(t)\)，等式两边同时取共轭转置有

\[ g^H(t) = y^H(t) P^H(t) = y^H(t) P(t) \tag{16} \]

此外由式（12）可导出增益向量的另一种计算式

\[ g(t) = h(t) / \alpha(t) \tag{17} \]

对上式两边同时取共轭转置有

\[ g^H(t) = h^H(t) / \alpha(t) = [h(t) / \alpha(t)]^H \tag{18} \]

将式（16）和式（18）带入式（15）有

\[ \begin{equation} \tag{19} \begin{aligned} W(t) &= W(t -1) - W(t -1) y(t) g^H(t) + x(t) g^H(t) \\ &= W(t -1) + [x(t) - W(t -1) y(t)] g^H(t) \\ &= W(t -1) + e(t) g^H(t) \\ &= W(t -1) + e(t) [h(t) / \alpha(t)]^H \end{aligned} \end{equation} \]

其中

\[ e(t) = x(t) - W(t -1) y(t) \tag{20} \]

整理以上各式得到 PAST 算法核心计算式如下：

\(y(i) = W^H(t) x(i)\)

\(h(t) = P(t - 1) y(t)\)

\(\alpha(t) = \beta + y^H(t) h(t)\)

\(g(t) = h(t) / \alpha(t)\)

\(P(t) = \beta^{-1} [P(t - 1) - h(t) h^H(t) / \alpha(t)]\)

\(e(t) = x(t) - W(t -1) y(t)\)

\(W(t) = W(t -1) + e(t) [h(t) / \alpha(t)]^H\)

上述过程也是 RLS 滤波器计算过程，而使用 RLS 方法来递归地计算 PAST，其关键一步是用 \(W^H(i -1) x(i)\) 来替代 \(W^H(t) x(i)\)，于是有 \(y(t) = W^H(t - 1) x(t)\)。

附 1：矩阵求逆引理

设 A 和 B 是两个 M x M 正定阵，它们之间的关系为

\[ \mathbf{A} = \mathbf{B}^{-1} + \mathbf{C} \mathbf{D}^{-1} \mathbf{C}^H \]

其中 D 是 N x N 正定阵，C 是 M x N 矩阵。根据矩阵求逆引理，A 的逆矩阵表示为

\[ \mathbf{A}^{-1} = \mathbf{B} - \mathbf{BC} \left(\mathbf{D} + \mathbf{C}^HBC \right)^{-1} \mathbf{C}^H \mathbf{B} \]

附 2：维纳滤波器

PAST 算法源于子空间分析中的极小化问题 \(\text{min} J(W)\) ，同时也可视为求解最优线性 FIR 滤波器问题，因此 PAST 算法的实现可利用维纳滤波器的结论。

维纳霍夫方程的矩阵形式表示为

\[ \mathbf{R w_o} = \mathbf{p} \]

其中 \(\mathbf{R} = \mathbb{E} [\mathbf{uu^H}]\) 为滤波器抽头输入向量的自相关矩阵，\(\mathbf{p} = \mathbb{E} [\mathbf{u}(n) d^*(n)]\) 为滤波器输入向量 \(\mathbf{u}\) 与期望响应 \(d(n)\) 的互相关向量。

最优线性 FIR 滤波器抽头权向量为

\[ \mathbf{w_o} = \mathbf{R}^{-1} \mathbf{p} \]

References

张贤达. 矩阵分析与应用（第2版）. 清华大学出版社, 2013.