记一下机器学习笔记 最小均方（LMS）算法

大家好，欢迎来到IT知识分享网。

于是这章书主要是简单介绍LMS算法的原理，并介绍几个简单的最优化方法，然后用物理热力学原理描述LMS算法的学习过程（这个部分太过高深只好跳过）

LMS滤波结构

原理上跟感知机也差不多，也是对包含一组共$M$个元素的$x_1,x_2,...,x_M$的输入用一个线性组合器处理，也就是对其进行加权求和，得出结果$y$，与期望响应$d$相比较，获得误差信号$e$，并由此修正权值，如下图：

这里比感知机还要简单的，直接将局部诱导域$v$作为输出$y$。因此可以表述成：

$$y(i) = w_1(i)x_1(i)+w_2(i)x_2(i)+...+w_M(i)x_M(i)=\sum_{k=1}^M w_k(i)x_k(i)$$

或者写成向量的形式：

$$y(i)=\textbf{x}(i)^T \textbf{w}(i)$$

$\textbf{w}(i)$即权值向量 $[w_1(i),w_2(i),...,w_M(i)]^T$， $i$表示迭代次数。
误差信号为期望响应跟输出的差，即：

$$e(i) = d(i)-y(i)$$

无约束最优化问题

LMS算法的目标就是找到一组权值向量，使其输出响应跟期望响应最接近。

设立一个代价函数$E(\textbf{w})$，其对权值向量连续可微，用来描述输出响应跟期望响应的差距，也就是值越小越好。于是我们的目标就是酱紫：
找到一个最优的权值向量$\textbf{w}^*$，对于任何$\textbf{w}$都有：

$$E(\textbf{w}^*)\le E(\textbf{w})$$
这是一个无约束最优化问题。其解决的一个必要条件就是 $\nabla E(\textbf{w}) = \textbf{0}$。
也就是：

$$\left[ \frac{ \partial E}{\partial w_1} , \frac{ \partial E}{\partial w_2} ,...,\frac{ \partial E}{\partial w_M} \right]^T = \textbf{0}$$

一般的解决方法是从一个初始权值向量$\textbf{w}(0)$开始，不断迭代产生新的权值向量$\textbf{w}(i)$，对于每一个权值向量其代价函数都要小于上一个的代价函数，即$E(\textbf{w}(i)) \lt E(\textbf{w}(i-1))$，如此往复直到代价函数足够小为止。或者说在一个M维的空间里，从一个点出发，不停地往代价函数减小的方向走，直到走到最低点。

最速下降法

也就是反向传播算法梯度下降的基本原理，在每一个位置$\textbf{w}(i)$求出当前位置的代价函数的梯度$\textbf g (i)$，再沿着梯度的反方向（正方向使代价函数增加）移动一段距离成为$\textbf{w}(i+1)$，也就是每次都顺着坡最陡的方向往下走一步。
梯度即为代价函数对权值向量的每一个元素求偏导：

$$\textbf g = \nabla E(\textbf w) = \frac { \partial E}{ \partial \textbf w }$$
权值向量的修正为：

$$\textbf w(i+1) = \textbf w(i) - \eta \textbf g(i)$$
$\eta$为一个标量，称为步长或学习率参数，可以理解为沿着梯度方向走的一步的大小。

理论上来说学习率参数$\eta$在足够小的时候，才能完全保证权值向量的修正是让代价函数一步比一步小的。但是$\eta$太小又会导致收敛速度过慢。

定义代价函数：

$$E(\textbf w) = \frac 1 2 \sum^N_{i=1}(d_i - y_i)^2 = \frac 1 2 \sum^N_{i=1}(d_i - \textbf w^T \textbf x_i)^2$$
那么就有：

$$\textbf g = \partial (\frac 1 2 \sum^N_{i=1}(d_i - \textbf w^T \textbf x_i)^2)/\partial \textbf w = - \sum^N_{i=1} e_i\textbf x_i$$

其中$e_i = d_i - \textbf w^T \textbf x_i$即误差值。于是权值向量的修正为：

$$\textbf w(i+1) = \textbf w(i) - \eta \textbf g(i) = \textbf w(i) + \eta \sum^N_{i=1} e_i\textbf x_i$$
N为样本数量。

> X x1 x2 [1,] 5. 11. [2,] -6. 4. [3,] 7. 4. [4,] -5. 8. [5,] -5. 6. [6,] 1. 1. [7,] 14.089176 -7. [8,] 3. -1. [9,] 6. -4. [10,] -10. 2. ...

R中%*%为矩阵相乘符号，t()为矩阵转置。

# X为点的坐标数据集，d为各点的正确分类，即期望响应，值为-1和1。 W = c(0,0) #初始化权值向量 eta = 1e-6 #学习率参数 n = 50 #迭代次数 MSE = c() #初始化均方差数组 for(i in 1:n){ y = X %*% W e = d - y #计算分类误差 MSE[i] = mean(e2) #记录每一步的均方差 W = W + eta * t(X) %*% e # 修正权值 } plot(MSE,type='l',xlab='iteration') #绘制均方差变化曲线 y = sign(X %*% W) qplot(x1,x2,color=factor(y)) #绘制分类结果

接下来将学习率参数$\eta$改为一个较大的值7.5e-6：

可见权值向量的轨迹从平滑变成了抖动。

牛顿法

$$\Delta E(\textbf w) = E(\textbf w(i)+\Delta \textbf w) - E(\textbf w(i)) = \textbf g^T(i)\Delta \textbf w+\frac 1 2 \Delta \textbf w^T \textbf H(i) \Delta \textbf w$$
其中 $\textbf H$为Hessian矩阵：

$$\textbf H = \begin{bmatrix} \frac{\partial^2 E}{\partial w_1^2} & \frac{\partial^2 E}{\partial w_1\partial w_2 }& \cdots&\frac{\partial^2 E}{\partial w_1\partial w_M}\\ \frac{\partial^2 E}{\partial w_2\partial w_1 } & \frac{\partial^2 E}{\partial w_2^2}&\cdots& \frac{\partial^2 E}{\partial w_2\partial w_M } \\ \vdots&\vdots&\ddots&\vdots\\\frac{\partial^2 E}{\partial w_M \partial w_1} & \frac{\partial^2 E}{\partial w_M\partial w_2 }& \cdots&\frac{\partial^2 E}{\partial w_M^2}\end{bmatrix}\quad$$
说白了就是对不同组合的权值求两次偏导。

接着就是要最大化$\Delta E(\textbf w)$，所以拿上上式右边对权值向量求导后再使之为0：

$$\textbf g(i)+\textbf H(i) \Delta \textbf w = \textbf 0$$
解得 $\Delta \textbf w = - \textbf H^{-1}(i)\textbf g(i)$。
也就是 $\textbf w(i+1) = \textbf w(i)- \textbf H^{-1}(i)\textbf g(i)$。

牛顿法的主要问题就是计算复杂度，以及其要求Hessian矩阵$\textbf H$每次迭代里都必须是正定的但这不好保证。

对于代价函数是这样的情况：

$$E(\textbf w) = \frac 1 2 \sum^N_{i=1}(d_i - y_i)^2 = \frac 1 2 \sum^N_{i=1}(d_i - \textbf w^T \textbf x_i)^2$$
拿代价函数对权值求两次偏导，可以算得Hessian矩阵 $\textbf H$的第i行第j列的元素为：

$$h_{ij}=-\sum^N_{s=1}x_i(s)x_j(s)$$
其中N为样本数量，s表示第s个样本。
因而Hessian就为：

$$\textbf H = \textbf X^T \textbf X$$
其中 $\textbf X$为样本矩阵，一行一样本一列一属性。

那么训练的R脚本就是酱紫：

H = t(X) %*% X #计算Hessian矩阵 W = c(0,0) #初始化权值向量 n=50 for(i in 1:n){ y = X %*% W e = d - y g = - t(X) %*% e W = W - solve(H) %*% g #按照公式修正权值 }

R中函数solve()可以求解矩阵的逆。
结果发现一次迭代就直接走到了最优值。

高斯-牛顿法

$$E(\textbf w) = \frac 1 2 \sum^N_{i=1}(d_i - y_i)^2 = \frac 1 2 \sum^N_{i=1}e_i^2$$
不过这次就不先拿代价函数，而是拿误差信号 $e_i$对权值向量在某一点处作一阶泰勒展开：

$$e'_i(\textbf w) = e_i(\textbf w(n))+\left[\frac{ \partial e_i}{ \partial \textbf w}\right]^T(\textbf w-\textbf w(n))$$
这回i表示第i个样本，而n表示第n次迭代。
把所有样本的 $e'_i$组合成列阵形式，就有：

$$\textbf e'(\textbf w) = \textbf e(\textbf w(n)) +\textbf J(n)(\textbf w-\textbf w(n))$$

其中$\textbf e' = \left[e_1',e_2',...,e_N'\right]^T$，N为样本数量。
$\textbf J$为Jacobi矩阵：

$$\textbf J = \begin{bmatrix} \frac{\partial e_1}{\partial w_1} & \frac{\partial e_1}{\partial w_2}& \cdots&\frac{\partial e_1}{\partial w_M}\\ \frac{\partial e_2}{\partial w_1}&\frac{\partial e_2}{\partial w_2}&\cdots& \frac{\partial e_2}{\partial w_M} \\ \vdots&\vdots&\ddots&\vdots\\ \frac{\partial e_N}{\partial w_1} & \frac{\partial e_N}{\partial w_2}& \cdots&\frac{\partial e_N}{\partial w_M}\end{bmatrix}\quad$$
说白了就是每个样本的误差信号分别对每个权值求偏导。

那么误差信号就是：

$$E(\textbf w) = \frac 1 2 \lVert\textbf e'(\textbf w)\rVert^2=\frac 1 2 \lVert\textbf e(\textbf w(n))\rVert^2 + \textbf e(\textbf w(n))^T \textbf J(n)(\textbf w-\textbf w(n)) + \frac 1 2(\textbf w-\textbf w(n))^T \textbf J(n)^T \textbf J(n)(\textbf w-\textbf w(n))$$
矩阵形式的完全平方公式。两根竖线 $\lVert$表示欧几里得范数，也就是常说的向量的模。
现在需要找到一个权值向量使上式最小作为 $\textbf w(n+1)$，于是对上式对权值向量求导并使之为0，得：

$$\textbf J^T(n)\textbf e(n)+\textbf J^T(n)\textbf J(n)(\textbf w-\textbf w(n))=\textbf 0$$
解出 $\textbf w$作为 $\textbf w(n+1)$得：

$$\textbf w(n+1) = \textbf w(n) -(\textbf J^T(n)\textbf J(n))^{-1}\textbf J^T\textbf e(n)$$
这就是高斯-牛顿法的基本型。

自然这里还有要求$\textbf J^T(n)\textbf J(n)$得是非负定的。于是通常会给它加上一个对角矩阵$\delta \textbf I$。$\delta$是一个较小的正数，$\textbf I$是单位矩阵。于是上式就变成：

$$\textbf w(n+1) = \textbf w(n) -(\textbf J^T(n)\textbf J(n)+\delta \textbf I)^{-1}\textbf J^T\textbf e(n)$$

维纳滤波器

然后接着推导。在这里误差信号为$e_i = d_i-y_i=d_i-\textbf w^T \textbf x_i$
于是有$\frac{ \partial e_i}{ \partial \textbf w} =- \textbf x_i$，$e'_i(\textbf w) = e_i(\textbf w(n))-\textbf x_i$。
从而Jacobi矩阵为：

$$\textbf J =\left[ - \textbf x_1,-\textbf x_2,...,-\textbf x_N \right]^T = -\textbf X$$
$\textbf X$就是样本矩阵。
另外可知有 $\textbf e = \textbf d - \textbf X \textbf w$。
将这些带入到高斯-牛顿法的基本型中可得：

$$\textbf w(n+1) = \textbf w(n) +(\textbf X^T(n)\textbf X(n))^{-1}\textbf X^T(\textbf d(n)- \textbf X \textbf w(n))$$
整理之后你会发现 $\textbf w(n)$会被消掉，然后就干脆成了：

$$\textbf w(n+1) = (\textbf X^T(n)\textbf X(n))^{-1}\textbf X^T(n)\textbf d(n)$$
简直可以一开始就一次计算啊。也难怪前边用牛顿法可以一次就收敛。
然后定义 $\textbf X$的伪逆为 $\textbf X^+=(\textbf X^T(n)\textbf X(n))^{-1}\textbf X^T$。这样就可以表述成最优权值向量为样本矩阵的伪逆乘上期望响应：

$$\textbf w(n+1) = \textbf X^+\textbf d(n)$$
这就像是《神机》第二章所讲的一次性计算分界的线性最小二乘分类器，所以这也叫 线性最小二乘滤波器。
当样本数量N趋于无穷时，就成了维纳滤波器。

R语言中用行代码即可算得权值向量：

W = solve(t(X) %*% X) %*% t(X) %*% d

最小均方算法

因此代价函数就成了：

$$E(\textbf w) = \frac 1 2 (d_i - y_i)^2 = \frac 1 2 e_i^2$$
称为代价函数的瞬时值。
求偏导后即可得：

$$\frac {\partial E}{\partial \textbf w(n)}=-\textbf x(n)\textbf e(n)$$
于是权值调整方式就是这样：

$$\textbf w(n+1)=\textbf w(n) + \eta\textbf x(n)\textbf e(n)$$
$\eta$同上为学习率参数。

R代码实现如下：

W = c(0,0) #初始化权值 eta = 1e-4 #学习率参数 n = 5 #进行5轮迭代 MSE = c() #均方根误差记录向量 for(t in 1:n){ for(i in 1:N){ e = d[i]-X[i,] %*% W W = W + e*X[i,] #修正权值 E = d - X %*% W MSE = append(MSE,mean(E2)) #计算并记录均方差 } }

这里学习率参数设为1e-4。可见在第一轮迭代中就已经收敛。因而在大量样本的数据中LMS的随机梯度方法相比前面几个方法更有性能优势。

这里是权值向量调整的轨迹。尽管是边抖边走也最终还是走到了最优处，到了目的地之后就在原地做起了布朗运动。

学习率退火方案

限制LMS算法性能的一个因素就是学习率参数$\eta$被设为是固定的，更科学的方式应该是一开始大，后面越来越小。
于是就有人提出了一个形式，学习率参数应该随迭代次数变化： $\eta(n)=\frac c n$。这里c是一个常数。
但是要是c设得比较大，导致一开始的时候$\eta$太大咋办？于是就又提出了下面的方式：

$$\eta(n)=\frac {\eta_0}{1+n/\tau}$$
这里 $\eta_0$和 $\tau$都是可调常数。酱紫就可以在一开始的时候 $\eta$不至于过大，而到了后期的时候接近于 $\frac c n$。这里 $c=\eta_0\tau$。