神经网络与机器学习
第5章 随机梯度下降法-BP的起源
神经网络的训练有很多方法，以数值优化为基础的随机梯度学习算法能够处理大规模的数据集合，它也是后面多层神经网络后向传播算法的基础。
随机梯度下降是以均方误差为目标函数的近似最速下降算法，该算法被广泛用于自适应信号处理领域，是Widrow与学生Hoff提出的，他们当时叫做自适应线性神经元（线性组合滤波器）和最小均方学习。
§5.1 自适应滤波器
我们先学习一下什么是滤波。滤波、平滑、预测都属于信号估计，三种方式的区别如下面图所示。滤波是给出信号$s(K)$在当前时刻$n$的估计$\widehat{s}(k)=f[x(k),\cdots,x(0)]$，$k$时刻前的数据可以用。
平滑是0在$k$到时刻数据中，$k$时刻数据缺少了，那么根据$k^'$前后时刻的数据估计$k^'$时刻的数据，称为平滑$\widehat{s}(k)=f[x(k),\cdots,x(k^'+1),x(k^'-1),\cdots,x(0)]$。数学里叫插值。
预测是依据从0到$k$时刻数据，得到并没有观测到$k+\tau$时刻的数据信息，给出一个估计$\widehat{s}(k+\tau)=f[x(k),\cdots,x(0)]$。
线性滤波器：是构造输入数据的函数$f[x(k),\cdots,x(0)]$，如果$f$为线性函数，比如
\[\widehat{s}(k)=\frac{x(k)}{4}+\frac{x(k-1)}{2}+\frac{x(k-2)}{4}\]
这个滤波器的结构非常简单，就是当前值为3个时刻观测值的加权平均值，线性系数[1/4，1/2，1/4]称为滤波器系数，这类线性函数是线性滤波器，否则是非线性滤波器$f$。维纳滤波和卡尔曼滤波都是线性滤波，维纳滤波适用于平稳的输入，而卡尔曼滤波是时变的最优线性滤波，适用于非稳态的输入。
自适应滤波：按照递归算法自动调节滤波器系数，平稳输入情况下收敛到维纳滤波器，非平稳情况下，如果变化足够缓慢，也能慢慢跟踪输入信号的变化，自适应滤波在通讯、地震、雷达、医学等方面应用广泛，主要可以完成系统辨识，均衡、预测、噪声对消这四个信号处理任务，应用十分广泛，并且由此衍生的随机梯度下降法对于机器学习也影响很大。
自适应滤波器是一类特殊的感知神经元，因为自适应滤波器由抽头延迟线、线性神经元组成的，如
从图5.1可以看出自适应横向滤波器就是一种特殊的线性神经元（感知机），激活函数选择的是线性函数$y=x$，我们可以将神经元$k$时刻的输出写为
\[a(k)=w^{\mathrm{T}}p+b=\sum_{i=1}^Rw_ip_i+b=\sum_{i=1}^Rw_ix(k-i+1)+b\]
如果增广向量
\[p\rightarrow \begin{bmatrix} p\\ 1 \end{bmatrix},w\rightarrow \begin{bmatrix} w\\ b \end{bmatrix}\] 5.3
自适应滤波器或者线性神经元可以简化为
$a(k)=\bm{w}^{\mathrm{T}}\bm{p}$ 5.4
自适应滤波器也和神经网络很像，有训练和测试两个过程，在训练过程中，已知期望响应$t(k)$，通过对权向量的设计，使从神经元输出在某种意义下最佳逼近期望响应。在工作过程中，权向量由训练过程得到，对未知输入（通常是随机信号）滤波得到响应信号。
性能曲面函数
误差信号
$e(k)=t(k)-a(k)=t(k)-\bm{W}^{\mathrm{T}}\bm{p}$ 5.5
其平方称为瞬时平方误差
$e^2(k)=t^2(k)-2t(k)\bm{w}^{\mathrm{T}}\bm{p}+\bm{w}^{\mathrm{T}}\bm{p}\bm{p}^{\mathrm{T}}\bm{w}$ 5.6
设$p(k)$、$t(k)$是平稳的随机过程样本-随机信号，满足一定的分布，那么我们希望均方误差
$E[e^2]=E[t^2]-2\bm{w}^{\mathrm{T}}E[t\bm{p}]+\bm{w}^{\mathrm{T}}E[\bm{p}\bm{p}^{\mathrm{T}}]\bm{w}$ 5.7
最小，这里输入数据向量$p$的自相关矩阵
\[R=E[pp^{\mathrm{T}}]\] 5.8
输入和期望信号的互相关向量为
$d=E[t\bm{p}]$ 5.9
那么可以进一步将均方误差$E[\epsilon^2(k)]$写成如下形式
\[\xi =E[e^2(k)]=E[t^2(k)]-2dw+w^{\mathrm{T}}Rw\] 5.10
可以看出均方误差是权向量$w$的二次函数，其二阶导数-Hessian矩阵
\[\triangledown ^2\xi =\frac{\partial ^2\xi}{\partial w^2} =2R\succeq 0\]
$R\succeq 0$表示$R$矩阵非负定，我们一般遇到$R\succ 0$正定，因此均方误差是权向量$w$的二次凸函数，必定有最小解，那么通过一次导数等于零，解出最优权向量$w^o$
\[\triangledown \xi =\frac{\partial \xi}{\partial w} =2Rw- 2d=0\Rightarrow w^o=R^{-1}d\] 5.11
那么最优权向量$w^o$代入公式（3），最小均方误差$\xi_{min}$写为
\[\xi_{min}=E[t^2(k)]-d^{\mathrm{T}}w^o=E[t^2(k)]-d^{\mathrm{T}}R^{-1}d\] 5.12
我们也一般利用第4章学的最速下降法、牛顿法等迭代求解最优权向量$w^o$。
§5.2 随机梯度下降法
我们现实中遇到这样的问题，自相关矩阵$R$得不到，我们不知道随机输入信号的分布，我们只知道样本集合
\[\left \{ p_1,t_1 \right \},\left \{ p_2,t_2 \right \},\cdots,\left \{ p_Q,t_Q \right \}\]
下标$q=1,2,\cdots,Q$是样本数量。最速下降法、牛顿法等迭代求解方法都无法进行，Widrow提出了一种近似最速下降算法，将均方误差$\xi=E[e^2(k)]$被瞬时平方误差代替
\[\widetilde{\xi}=e^2(k)\]
那么每次迭代计算寻优时，都计算梯度估计值
\[\widetilde{\triangledown\xi}=\triangledown e^2(k)=2e(k)\triangledown e(k)=2e(k)\frac{\partial e(k)}{\partial w}\]
由于$e(k)=t(k)-a(k)=t(k)-\bm{W}^{\mathrm{T}}\bm{p}$，所以上式可以写为
$\widetilde{\triangledown\xi}=\triangledown e^2(k)=-2e(k)\bm{p}$
可以看出梯度近似值只是误差和输入向量的乘积而已，所以依次放入样本${ p_1,t_1},{ p_2,t_2 },\cdots,{ p_Q,t_Q }$，恒定的学习率$\alpha$，则依次更新权向量
\[w(q+1)=w(q)-\alpha \widetilde{\triangledown\xi}=w(q)+2\alpha e(k)p\]
这样每次样本输入，神经网络权值进行更新，这种学习方式又称为在线学习，最小均方算法，$delta$规则，随机梯度下降法，或者Widrow-Hoff学习。
按照第4章最优化中的二次函数分析，理论分析学习率$\alpha<\frac{2}{\lambda_{max}}$与自相关矩阵$R$最大特征值$\lambda_{max}$相关，但是现实中不知道自相关矩阵，因此学习率经常是试错法选择。理论分析学习可以在自适应滤波理论中详细学习，这里不再一一讲解。
这里列举两个例子来看随机梯度下降法的应用。一是利用随机梯度下降法训练感知机，另外一个是自适应噪声对消。
例子5.1 感知机的训练
我们有两个样本分别表示橘子和苹果，各自的期望目标输出为
\[\left \{ p_1=\begin{bmatrix} 1\\ -1\\ -1 \end{bmatrix},t_1=-1 \right \},\left \{ p_2=\begin{bmatrix} 1\\ -1\\ -1 \end{bmatrix},t_2=-1 \right \}\]
我们利用随机梯度下降法，依次输入${p_1,t_1},{p_2,t_2}$，${p_1,t_1},{p_2,t_2}$，...,因为只有两个样本，所以可以组成这样的在线训练序列，${p_1,t_1},{p_2,t_2}$的顺序可以随机输入。初始化，我们选择这样的感知机
\[a=w^{\mathrm{T}}p=[w_1,w_2,w_3]\begin{bmatrix} p_1\\ p_2\\ \vdots\\ p_R \end{bmatrix}\]
$f$取线性函数$f(x)=x$ ，偏置为$b=0$，初始的权向量为$w(0)=[0,0,0]^{\mathrm{T}}$，学习率$\alpha=0.2$，那么我们输入${ p_1=\begin{bmatrix} 1\\ -1\\ -1 \end{bmatrix},t_1=-1 }$，得到
\[a(0)=w^{\mathrm{T}}(0)p(0)=w^{\mathrm{T}}(0)p_1=[0,0,0]\begin{bmatrix} 1\\ -1\\ -1 \end{bmatrix}=0\]
误差为
\[e(0)=t(0)-a(0)=t_1-a(0)=-1\]
那么按照随机梯度下降方法更新权
\[w(1)=w(0)+2\alpha e(0)p_1=\begin{bmatrix} 0\\ 0\\ 0 \end{bmatrix}+2\times 0.2\times (-1)\times \begin{bmatrix} 1\\ -1\\ -1 \end{bmatrix}=\begin{bmatrix} -0.4\\ 0.4\\ 0.4 \end{bmatrix}\]
我们输入${ p_2=\begin{bmatrix} 1\\ -1\\ -1 \end{bmatrix},t_2=-1 }$，得到
\[a(1)=w^{\mathrm{T}}(1)p(1)=w^{\mathrm{T}}(1)p_2=[-0.4,0.4,0.4]\begin{bmatrix} 1\\ 1\\ -1 \end{bmatrix}=-0.4\]
误差为
\[e(1)=t(1)-a(1)=t_2-a(0)=1.4\]
那么按照随机梯度下降方法更新权
\[w(2)=w(1)+2\alpha e(1)p_2=\begin{bmatrix} -0.4\\ 0.4\\ 0.4 \end{bmatrix}+2\times 0.2\times (1.4)\times \begin{bmatrix} 1\\ 1\\ -1 \end{bmatrix}=\begin{bmatrix} 0.16\\ 0.96\\ -0.16 \end{bmatrix}\]
我们再次输入${p_1,t_1}$，得到
\[a(2)=w^{\mathrm{T}}(2)p(2)=w^{\mathrm{T}}(2)p_1=-0.64\]
\[e(2)=t(2)-a(2)=t_1-a(2)=-0.36\]
\[w(3)=w(2)+2\alpha e(2)p_1=\begin{bmatrix} 0.016\\ 1.104\\ -0.016 \end{bmatrix}\]
这样训练下去，我们会得到
\[w(\infty)=\begin{bmatrix} 0\\ 1\\ 0 \end{bmatrix}\]
这比前面第3章得到感知机决策边界更加好，因为前面的决策边界模式一旦划分正确则停止，而现在的随机梯度下降法是使得决策边界尽量原来样本模式，使得近似均方误差最小，两种决策边界虽然都处于样本模式之间，但是随机梯度法得到决策边界在测试中是非常实用的。
例子5.2 自适应噪声对消
图5.3 自适应噪声对消
任务：设噪声源$u(k)$将语音等信号$s(k)$污染，我们不知道噪声强度多少，我们想在语音源附近也测量同源噪声$\beta(k)$，然后用自适应算法预测噪声$u(k)$，然后对消噪声，得到误差就是相对纯净的信号发送出去。
设信号是随机相位信号
\[s(k)=cos(\frac{\pi k}{16}+\varphi )\]
污染信号的噪声是
\[u(k)=0.9u(k)+\eta(k)\]
$\eta(k)$是标准正态分布的高斯白噪,$\sigma_{\eta}^2=1$，而另外一个渠道得到的噪声和$u(k)$同源，但是不同的信道
\[\beta(k)=-0.75\beta(k)+\eta(k)\]
用随机梯度下降算法进行自适应滤波，不需要任何关于随机输入的二阶矩的知识。主程序如下：
%System parameters
omg0 = pi/16; a1 = 0.9; a2 = -0.75; var_eta = 1;
% Generate Desired and Input signals
order=20;
N = 1000;%data length
n = [0:N-1]';
phi = 2*pi*rand(1,1);
s = cos(omg0*n+phi);
eta = randn(N,1);
u = filter(1,[1,-a1],eta);
beta = filter(1,[1,-a2],eta);
t = s + u;
mu = 0.005;%step
w0 = zeros(order,1);
[w,e,yhat] = firlms(beta,t,mu,order,w0);
plot(e,'b-');
hold on
plot(s,'r--')
| function [w,e,yhat] = firlms(x,t,mu,L,w0) % Computes filter coefficients w and the corresponding error in % e, given signal x, desired signal y, step size mu, filter order L, and % the initial coefficient vector w0. N = length(x); % Length of x(n) & number of iterations x = reshape(x,N,1); % x as a column vector w = zeros(L,N); % Initialization of w X = zeros(L,N); % Data matrix xx = [zeros(L-1,1);x]; % Prepend L-1 zeros for i = 1:L X(i,:) = xx(L-i+1:end-i+1)'; end e = zeros(1,N); yhat = zeros(N,1); %--LMS Algorithm for FIR filter - real-valued case % Algorithm initialization yhat(1) = w0'*X(:,1); e(1) = t(1) - yhat(1); w(:,1) = w0 + 2*mu*X(:,1)*e(1); % Iterations for n = 2:N yhat(n) = w(:,n-1)'*X(:,n); e(n) = t(n) - yhat(n); w(:,n) = w(:,n-1) + 2*mu*X(:,n)*e(n); end |
一开始滤波器并不能消除噪声，随着时间推移，慢慢滤波器系数（神经元权向量收敛），对消了噪声，余下相对纯净的信号。这种系统在飞机、战场、大风等恶劣天气中个人通信等等场景中经常遇到。尽管飞行器发动机、风啸等噪声和人的语音同时进入声音发射器（麦克风），但是相对清晰的语音被发送过去，当然如果遇到与语音频率特别接近的声音是无法消除的。
图5.4 自适应噪声对消示例
例子5.3（*选讲）：二阶自回归AR(2)过程
\[x(n)-0.96x(n_1)+0.9x(n-2)=\eta(n)\]
$\eta(n)$是高斯白噪声，均值为0，方差为$\sigma_{\eta}^2$，我们利用一个二阶滤波器进行系统辨识
\[y(n)=w_1x(n-1)+w_2x(n-2)\]
%AR过程
N = 501; mu= 0.01;
a1 = -0.950; a2 = 0.9; varx = 1; r0 = varx;
vareta = ((1-a2)*((1+a2)^2-a1^2)/(1+a2))*varx;
r1 = -a1/(1+a2)*r0; r2 = (-a2+a1^2/(1+a2))*r0;
lam1 = (1-a1/(1+a2))*varx; lam2 = (1+a1/(1+a2))*varx;
XR = lam1/lam2;d = [r1;r2];%条件数
R = toeplitz([r0,r1,r2]);%自相关阵
w = sqrt(vareta)*randn(N,1);
x = filter(1,[1 a1 a2],w);%传递函数形式
%观测随机梯度算法的性能
[w,e] = lplms(x,0.01,2,[0 0]');
subplot(121)
plot(w(1,:),w(2,:));
hold on
[w1,w2]=meshgrid(-1:0.01:3,-2:0.01:1);
z=r0+r0*w1.^2+r0*w2.^2+2*r1*w1.*w2-2*r1*w1-2*r2*w2;%MSE
contour(w1,w2,z,50);
subplot(122)
plot(w(1,:),'b--*');
hold on
plot(w(2,:),'r-^')
plot([0 501],[0.95 0.95],'--')
plot([0 501],[-0.9 -0.9],'--')
| function [w,e] = lplms(x,mu,L,w0) % Computes filter coefficients in w and the corresponding error in % e, given signal x, desired signal d, step size mu, filter order L, and % the initialcoefficient vector w0. N = length(x); % Length of x(n) x = reshape(x,N,1); % x as a column vector d = x(2:end); % Generate desired d(n) N = length(d); % Number of iterations c = zeros(L,N); % Initialization of c X = zeros(L,N); % Data matrix xx = [zeros(L-1,1);x]; for i = 1:L X(i,:) = xx(L-i+1:end-i)'; end e = zeros(1,N); %--LMS Algorithm % Algorithm initialization yhat = w0'*X(:,1); e(1) = d(1) - yhat; w(:,1) = 2*mu*X(:,1)*e(1); % Iterations for n = 2:N yhat = w(:,n-1)'*X(:,n); e(n) = d(n) - yhat; w(:,n) = w(:,n-1) + 2*mu*X(:,n)*e(n); end w = [w0,w]; e = [x(1),e]; |
图5.5 随机梯度下降法中函数值和权向量收敛学习曲线
\[y(n)=w_1x(n-1)+w_2x(n-2)\]
权系数$w_1$和$w_2$相对"随机"地收敛到-0.95和0.9，达到系统辨识的目的。
作业：利用随机梯度下降算法训练四分类判定决策界面问题。
图3.2 四类数据的几何表示
训练数据为8个目标向量，分别是
\[\left \{ (p_1=\begin{bmatrix} 1\\ 1 \end{bmatrix},t_1=\begin{bmatrix} -1\\ -1 \end{bmatrix}),(p_2=\begin{bmatrix} 1\\ 2 \end{bmatrix},t_2=\begin{bmatrix} -1\\ -1 \end{bmatrix}) \right \}\]
类1标记为○
\[\left \{ (p_3=\begin{bmatrix} 2\\ -1 \end{bmatrix},t_3=\begin{bmatrix} -1\\ 1 \end{bmatrix}),(p_4=\begin{bmatrix} 2\\ 0 \end{bmatrix},t_4=\begin{bmatrix} -1\\ 1 \end{bmatrix}) \right \}\]
类2标记为□
\[\left \{ (p_5=\begin{bmatrix} -1\\ 2 \end{bmatrix},t_5=\begin{bmatrix} 1\\ -1 \end{bmatrix}),(p_6=\begin{bmatrix} -2\\ 1 \end{bmatrix},t_6=\begin{bmatrix} 1\\ -1 \end{bmatrix}) \right \}\]
类3标记为●
\[\left \{ (p_7=\begin{bmatrix} -1\\ -1 \end{bmatrix},t_7=\begin{bmatrix} 1\\ 1 \end{bmatrix}),(p_8=\begin{bmatrix} -2\\ -2 \end{bmatrix},t_8=\begin{bmatrix} 1\\ 1 \end{bmatrix}) \right \}\]
类4标记为■
神经元选为带偏置的
\[a=Wp+b=\begin{bmatrix} W & b \end{bmatrix}\begin{bmatrix} p\\ 1 \end{bmatrix}\]
设初值
\[W(0)=\begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix},b(0)=\begin{bmatrix} 1\\ 1 \end{bmatrix}\]
或者增广
\[W(0)=\begin{bmatrix} 1 & 0 & 1\\ 0 & 1 & 1 \end{bmatrix}\]
那么利用随机梯度下降算法40次，并画出收敛的决策边界。
| 2021年3月27日星期六 clc,clear |