学习篇：自适应算法总结（LS、LMS、TLS）

本内容为博主在学习TLS算法的过程中的一些总结对比，资料来源于张贤达的《矩阵分析》以及《现代信号处理》，部分来源于网络。

自适应滤波首先得从维纳滤波讲起，LS以及LMS本质上都是因为不知道统计特性用一段时间内的数据来代替期望值，本质还是要收敛到维纳解。

1. 维纳滤波

因果维纳滤波器的输出 $y(n)$ :

y(n) = \hat{s}(n) = x(x) * h(n) = \sum_{k=0}^{+\infin}{h_kx(n-k)}

设期望信号为 $d(n)$ ，误差信号 $e(n)$ ，以及均方值 $E\left[|e(n)|^2\right]$ 分别表示为：

e(n)=d(n)-y(n)=s(n)-y(n)

代价函数为 $\Rightarrow J(n)=E\left[|e(n)|^2\right]=E\left[e(n) e^*(n)\right]$

设滤波器参数 $h(n)$ 为：

h(k)=a_k+j b_k, k=0,1,2, \ldots

求代价函数的偏导，使均方误差最小，满足：

\begin{aligned} & \min _{h_k} J(n) \rightarrow \nabla J(n)=\frac{\partial J(n)}{\partial h_k}=0 \\ & \frac{\partial E\left[|e(n)|^2\right]}{\partial a_k}+j \frac{\partial E\left[|e(n)|^2\right]}{\partial b_k}=0 \quad k=0,1,2, \ldots \end{aligned}

记梯度算子为

\begin{gathered} \nabla_k=\frac{\partial}{\partial a_k}+j \frac{\partial}{\partial b_k} \quad k=0,1,2, \ldots \\ \nabla_k J(n)=\frac{\partial E\left\{e(n) e^*(n)\right\}}{\partial a_k}+j \frac{\partial E\left\{e(n) e^*(n)\right\}}{\partial b_k} \end{gathered}

上式展开为

\nabla_k E\left[|e(n)|^2\right]=E\left[\frac{\partial e(n)}{\partial a_k} e^*(n)+\frac{\partial e^*(n)}{\partial a_k} e(n)+j \frac{\partial e(n)}{\partial b_k} e^*(n)+j \frac{\partial e^*(n)}{\partial b_k} e(n)\right]

又

\begin{aligned} e(n)=s(n)-y(n) & =s(n)-\sum_{k=0}^{+\infty} h_k x(n-k) \\ & =s(n)-\sum_{k=0}^{+\infty}(a(k)+j b(k)) x(n-k) \end{aligned}

\begin{aligned} & \frac{\partial e(n)}{\partial a_k}=-x(n-k) \\ & \frac{\partial e(n)}{\partial b_k}=-j x(n-k) \\ & \frac{\partial e^*(n)}{\partial a_k}=-x^*(n-k) \\ & \frac{\partial e^*(n)}{\partial b_k}=j x^*(n-k) \end{aligned}

将上述(8)式代入得

\nabla_k J(n)=\nabla_k E\left[|e(n)|^2\right]=-2 E\left[x^*(n-k) e(n)\right]

即：

E\left[x^*(n-k) e(n)\right] = 0

分析：如果要使滤波器的均方误差达到最小值，则误差信号和输入信号正交，也成为正交性原理。

由（10）式，可得：

\begin{aligned} & E\left[x^*(n-k) e_{o p t}(n)\right]=0, k=0,1,2, \ldots \\ & \Rightarrow E\left[x(n-k)\left(d(n)-\sum_{i=0}^{+\infty} h_{o p t, i} x(n-i)\right)^*\right\rceil=0 \end{aligned}

将输入信号分配进去, 得到

r_{x d}(-k)=\sum_{i=0}^{+\infty} h_{o p t, i}^* r_{x x}(i-k) \quad k=0,1,2, \ldots

维纳-霍夫 (Wiener-Hopf）方程:

r_{x d}(k)=\sum_{i=0}^{+\infty} h_{o p t, i} r_{x x}(k-i) \quad k=0,1,2, \ldots

列出N个时刻的等式，改写为矩阵形式得：

R_{xd} = R_{xx}h，即 h = h_{opt} = R^{-1}_{xx}R_{xd}

2. LMS自适应滤波

针对于上面的维纳解，我们需要知道 $R_{xx}，R_{xd}$ ，而且矩阵求逆计算量非常大，同时当信号非平稳时 $R_{xx}，R_{xd}$ 还可能随时间而变化，这种情况下我们可以用递归迭代的算法进行求解。

设误差信号为:

e(n)=d(n)-y(n)=d(n)-x^T(n) w(n)=d(n)-\boldsymbol{w}^T(n) \boldsymbol{x}(n)

均方误差为:

\begin{aligned} & J=E\left[e^2(n)\right]=E\left[d^2(n)\right]+\boldsymbol{w}^T(n) E\left[\boldsymbol{x}(n) \boldsymbol{x}^T(n)\right] \boldsymbol{w}(n) \\ & -2 E\left[d(n) \boldsymbol{x}^T(n)\right] \boldsymbol{w}(n) \end{aligned}

LMS自适应算法直接利用瞬态均方误差对瞬时抽头向量（滤波器系数) 求梯度:

\hat{\nabla} \boldsymbol{J}=\frac{\partial\left(e^2(n)\right)}{\partial w(n)}=\left[2 e(n) \frac{\partial e(n)}{\partial w_1(n)}, 2 e(n) \frac{\partial e(n)}{\partial w_2(n)}, \cdots, 2 e(n) \frac{\partial e(n)}{\partial w_M(n)}\right]^T=-2 e(n) \boldsymbol{x}(n)

由此可得传统LMS自适应滤波算法流程如下:

\begin{aligned} & y(n)=\boldsymbol{w}^T(n) \boldsymbol{x}(n) \\ & \boldsymbol{e}(n)=d(n)-y(n) \\ & \boldsymbol{w}(n+1)=\boldsymbol{w}(n)+2 \mu e(n) \boldsymbol{x}(n) \end{aligned}

3. LS自适应滤波（RLS）

LS一般都称作最小二乘法，自适应滤波时常采用RLS算法，即递归最小二乘法。（可以理解为LS解的一种迭代求法，否则和维纳解一样是需要矩阵求逆）

对比LMS的代价函数：

J_{LMS} = E\left\{|d(n)-\bf{w^Hu(n)}|^2\right\} \\ J_{LS} = \sum_i|d(i)-\bf{w^Hu(i)}|^2

也就是说，其实本质都是一样的，一个是求和，一个是求期望，如果 $i$ 的取值够大，那么两个应该是一样的。我们同样对上述的代价函数求梯度，得到：

\mathrm{J}=\|\overrightarrow{\mathrm{e}}\|^2=\vec{e}^{\mathrm{H}} \overrightarrow{\mathrm{e}}=(\overrightarrow{\mathrm{b}}-\mathrm{A} \overrightarrow{\mathrm{w}})^{\mathrm{H}}(\overrightarrow{\mathrm{b}}-\mathrm{A} \overrightarrow{\mathrm{w}})

\mathrm{J}=\overrightarrow{\mathrm{b}}^{\mathrm{H}} \overrightarrow{\mathrm{b}}-\overrightarrow{\mathrm{b}}^{\mathrm{H}} \mathrm{A} \overrightarrow{\mathrm{w}}-\overrightarrow{\mathrm{w}}^{\mathrm{H}} \mathrm{A}^{\mathrm{H}} \overrightarrow{\mathrm{b}}+\overrightarrow{\mathrm{w}}^{\mathrm{H}} \mathrm{A}^{\mathrm{H}} \mathrm{A} \overrightarrow{\mathrm{w}}

其中A是输入矩阵，b是期望值：

\nabla \mathrm{J}=-2 \mathrm{~A}^{\mathrm{H}} \overrightarrow{\mathrm{b}}+2 \mathrm{~A}^{\mathrm{H}} \mathrm{A} \overrightarrow{\mathrm{w}}

\mathrm{A}^{\mathrm{H}} \mathrm{A} \overrightarrow{\hat{w}}=\mathrm{A}^{\mathrm{H}} \overrightarrow{\mathrm{b}}

可以得到方程的解:

\overrightarrow{\hat{w}}=\left(\mathrm{A}^{\mathrm{H}} \mathrm{A}\right)^{-1} \mathrm{~A}^{\mathrm{H}} \overrightarrow{\mathrm{b}}

$\mathrm{A}^{\mathrm{H}} \mathrm{A}$ 就是维纳解中的 $R_{xx}$ , $A^H\overrightarrow{b}$ 就是 $R_{xd}$ 。

所谓的最小二乘估计，是指在每个时刻对所有已输入的信号而言，误差平方和最小。LS用于接收到的数据块长度一定，并且数据、噪声（干扰）的统计特性未知或者非平稳的情况，其优化目标是使得基于该数据块的估计与目标数据块间加权的欧几里德距离最小。

3. 总体最小二乘法（TLS）

TLS算法可以理解为LS算法的一种改进，在LS算法中，我们考虑的误差 $e = |d(i)-\bf{w^Hu(i)}|$ ，实际上的意思就是我确认我每个时刻的输入 $u(i)$ 是准确的，偏差只出现在拟合中；而TLS算法则是考虑输入输出都存在噪声的情况，即我观测到的输入也可能是不准确的。这两种差别可以通过下面的直线拟合图来体会：

其中，虚线表示LS的误差，实线表示的是TLS的误差。可以直观的感受到其实TLS的拟合效果会更好，使用的是点到直线的距离。实际上，这也比较符合我们实际应用，因为输入数据是我们测量得到的，那么就一定会存在测量误差，使用TLS算法的最优化问题就是最小化输入误差和输出误差的和。

问题提出：给定数据向量 $b$ 和一个数据矩阵 $A$ ，考虑一个线性方程 $Ax = b$ ，求 $x$ 的最优解。

我们在以前求解这个问题时，前提条件是 $b$ 收到了高斯白噪声干扰，采用最小二乘法能够使得误差平方和最小，称为最小方差无偏估计（MVU）。

\bf{\hat{x}}=\bf{\left(\mathrm{A}^{\mathrm{H}} \mathrm{A}\right)^{-1} \mathrm{~A}^{\mathrm{H}} \mathrm{b}}

如果说，A 存在误差或者扰动，那么最小二乘估计从统计观点看就不是最优的，它将是有偏的。

上面的求解等价于：

在条件\quad Ax = b + \Delta b \quad 约束下，\underset{x}{min}||\Delta b||

而TLS考虑的问题是：

\bf{(A + E)x = b + e} \\ \operatorname{TLS}: \quad \min _{\Delta \mathbf{A}, \Delta \mathbf{b}, \mathbf{x}}\|\Delta \mathbf{A}\|_{\mathrm{F}}^2+\|\Delta \mathbf{b}\|_2^2 \\ \text { subject to }(\mathbf{A}+\Delta \mathbf{A}) \mathbf{x}=\mathbf{b}+\Delta \mathbf{b}

其中 $\mathrm{F}$ 表示Frobineous范数。注意，约束条件 $(\mathbf{A}+\Delta \mathbf{A}) \mathbf{x}=\mathbf{b}+\Delta \mathbf{b}$ 有时也表示为 $(\mathbf{b}+\Delta \mathbf{b}) \in \operatorname{Range}(\mathbf{A}+\Delta \mathbf{A})$ 。

一种求解方法

根据形成的TLS问题，我们可以构建如下的拉格朗日函数:

\mathcal{L}(\Delta \mathbf{A}, \Delta \mathbf{b}, \boldsymbol{\lambda})=\|\Delta \mathbf{A}\|_{\mathrm{F}}^2+\|\Delta \mathbf{b}\|_2^2+2 \boldsymbol{\lambda}((\mathbf{A}+\Delta \mathbf{A}) \mathbf{x}-\mathbf{b}-\Delta \mathbf{b})

通过利用KKT条件，我们知道，当且仅当存在 $\lambda \in \mathbb{R}^m$ 使得

\begin{aligned} & \nabla \mathcal{L}_{\Delta \mathbf{A}}=2 \Delta \mathbf{A}+2 \boldsymbol{\lambda} \mathbf{x}^{\top}=\mathbf{0} \\ & \nabla \mathcal{L}_{\Delta \mathbf{b}}=2 \Delta \mathbf{b}-2 \boldsymbol{\lambda}=\mathbf{0} \\ & (\mathbf{A}+\Delta \mathbf{A}) \mathbf{x}=\mathbf{b}+\Delta \mathbf{b} \end{aligned}

结合上述条件(1)和(2)，可知 $\Delta \mathbf{A}=-\Delta \mathbf{b} \mathbf{x}^{\top}$ ，将该结果带入上式(3)可知，因此

\Delta \mathbf{b}=\frac{\mathbf{A x}-\mathbf{b}}{\|\mathbf{x}\|^2+1}

随后，将该结果带入 $\Delta \mathbf{A}=-\Delta \mathbf{b x}$ ，因此

\Delta \mathbf{A}=-\frac{(\mathbf{A x}-\mathbf{b}) \mathbf{x}^{\top}}{\|\mathbf{x}\|^2+1}

随后，我们将 $\Delta \mathbf{b}$ 和 $\Delta \mathbf{A}$ 带入目标函数，因此目标函数可以写作

\min _{\Delta \mathbf{A}, \Delta \mathbf{b}, \mathbf{x}} \frac{(\mathbf{A} \mathbf{x}-\mathbf{b})^2}{\|\mathbf{x}\|^2+1}

上面的推导过程也可以写成：

([-b, A]+[-e, E])\left[\begin{array}{l} 1 \\ x \end{array}\right]=0

或等价为

(B+D) z=0

式中, 增广矩阵 $\boldsymbol{B}=[-\boldsymbol{b}, \boldsymbol{A}]$ 和扰动矩阵 $\boldsymbol{D}=[-\boldsymbol{e}, \boldsymbol{E}]$ 均为 $m \times(n+1)$ 维矩阵, 而 $z=\left[\begin{array}{l}1 \\ x\end{array}\right]$ 为 $(n+1) \times 1$ 向量，有：

\bf{Bz = r =-Dz}

$r$ 可以被视为矩阵方程 $Bz = 0$ 的总体最小二乘解 z的误差向量。

约束条件为：

z^Hz = 1

注：此处明显z的模>1，也可以假设模等于任意大于1的常数，在下面的求偏导中都会消去，不影响结果。

优化问题可以写为：

\bf{min ||Bz||_2^2 = min ||r||^2}

定义目标函数为

J=\|\boldsymbol{B} \boldsymbol{z}\|_2^2+\lambda\left(1-\boldsymbol{z}^{\mathrm{H}} \boldsymbol{z}\right)

式中, $\lambda$ 为 Lagrange 乘数。注意到 $\|\boldsymbol{B} \boldsymbol{z}\|_2^2=\boldsymbol{z}^{\mathrm{H}} \boldsymbol{B}^{\mathrm{H}} \boldsymbol{B} \boldsymbol{z}$ , 故由 $\frac{\partial J}{\partial \boldsymbol{z}^*}=0$ , 得到

\boldsymbol{B}^{\mathrm{H}} \boldsymbol{B} \boldsymbol{z}=\lambda \boldsymbol{z}

这表明, Lagrange 乘数应该选择为矩阵 $\boldsymbol{B}^{\mathrm{H}} \boldsymbol{B}$ 的最小特征值 (即 $\boldsymbol{B}$ 的最小奇异值的平方桹), 而总体最小二乘解 $\boldsymbol{z}$ 是与最小奇异值 $\sqrt{\lambda}$ 对应的右奇异向量。
令 $m \times(n+1)$ 增广矩阵 $\boldsymbol{B}$ 的奇异值分解为