《神经网络与机器学习》第7讲后向传播算法局限性与变种算法
神经网络与机器学习
第7章 反向传播算法局限性和变形
为什么反向传播算法可以训练神经网络，这是因为假如权值$w_{i,j}^m$做了一个小小的改动$\Delta w_{i,j}^m$改变，会引起性能函数的改变$\Delta F$
图7.1 多层前馈神经网络
\[\Delta F=\frac{\partial F}{\partial w_{i,j}^m}\Delta w_{i,j}^m\]
这个改变的比例就是我们求的链式法则的导数部分，因此反向传播就是细致地追踪一个$w_{i,j}^m$的微小变化如何导致性能函数的改变，反之性能函数的改变也反之作用到权值的改变上。
但是反向传播算法具有很多局限性，第一个是收敛性问题，其次有梯度消失和爆炸问题，还有过拟合，泛化能力差等问题。这一章我们主要探讨一下收敛性和改进的变形算法。
§7.1 非凸的前馈多层神经网络性能曲面
比如我们考虑下面一个$1\times 2\times 1$网络
图7.2 网络
我们已知最优解向量
\[w^1=\begin{bmatrix} w_{1,1}^1\\ w_{2,1}^1 \end{bmatrix}=\begin{bmatrix} 10\\ 10 \end{bmatrix},b^1=\begin{bmatrix} b_1^1\\ b_2^1 \end{bmatrix}=\begin{bmatrix} -5\\ 5 \end{bmatrix},w^2=\begin{bmatrix} w_{1,1}^2 & w_{1,2}^2 \end{bmatrix}=\begin{bmatrix} 1 & 1 \end{bmatrix},b^2=-1\]
两层的神经元函数都为logsig函数
\[f^1(x)=f^2(x)=\frac{1}{1+e^{-x}}\]
输入区间为$p\in [-2,2]$，设此时的函数值是真实的，虽然人为设计的一个函数，但是说明下面一些问题。
图7.3 目标函数
图7.3 对应权系数$w_{1,1}^1$和$w_{1,1}^2$的平方误差曲面和等高线
利用下面这段程序可以看到对应权系数$w_{1,1}^1$和$w_{1,1}^2$的平方误差曲面和等高线，见图7.3。它明显不是二次曲面，曲率在空间中变化多端，局部最小点各种各样，所以基于最速下降的随机梯度算法难免会陷入某个局部最小点。
%设变量为 W1_{1,1}, W2_{1,1}
[X,Y]=meshgrid(-5:0.5:15);
Z=arrayfun(@(x,y) fcurve(x,y),X,Y);
surf(X,Y,Z)
figure
contour(X,Y,Z,50)
function y=fcurve(w1,w2)
W1 =[10;10];
b1 = [-5;5];
W2 = [1 1];
b2 = -1;
P=[-2:0.1:2];
for k=1:length(P)
p=P(k);%
a1 = 1./(1+exp(-(W1*p+b1)));
a2 = 1./(1+exp(-(W2*a1+b2)));
T(k)=a2;
a1 = 1./(1+exp(-([w1;10]*p+b1)));
a2 = 1./(1+exp(-([w2 1]*a1+b2)));
Tw(k)=a2;
end
y=sum((T-Tw).^2)/length(T);
end
同样地，如果仅仅调节$w_{1,1}^1$和$b_1^1$的平方误差曲面，如图7.4。误差曲面更加扭曲，有的地方非常陡峭，而其他区域则非常平缓，随机梯度下降算法则会遇到困难，它在平坦的地方则停止不前，不会收敛到最小值处，这带来很大的挑战。
图7.4 对应权系数$w_{1,1}^1$和偏置$b_1^1$的平方误差曲面和等高线
同样地，如果仅仅调节$b_1^1$和$b_2^1$的平方误差曲面，如图7.5。误差曲面是对称的，这里出现鞍点。这里提示我们经常初始参数不取零点，因为可能是一个鞍点。所以初始参数经常取随机值，但是不能太大，因为会跑到其他局部最小值或者很平缓的地方。
图7.5 对应偏置系数$b_1^1$和$b_2^1$的平方误差曲面和等高线
§7.2 批处理batching的随机梯度下降算法
第6章中随机梯度下降算法是在线算法，每一个样本经过网络，则网络相应地改变一次，这里我们介绍批处理训练方法，设一组样本含有$Q$个，设每个样本等概率出现，那么网络的均方误差性能指标可以写为
\[F=\frac{1}{Q}\sum_{q=1}^Q(t_q-a_q)^{\mathrm{T}}(t_q-a_q)\]
注意每个样本都对应着网络输出$a_q$和目标值$t_q$都是向量，那么总的梯度为
\[\triangledown F=\frac{1}{Q}\sum_{q=1}^Q\triangledown (t_q-a_q)^{\mathrm{T}}(t_q-a_q)\]
就是各个样本的梯度平均值。那么我们把所有样本的前向计算和敏感度都计算出来，得到各个权系数的梯度平均，那么批处理的随机梯度下降算法更新公式
\[\left\{\begin{matrix} W^m(k+1)=W^m(k)-\frac{\alpha}{Q}\sum_{q=1}^Q\delta_q^m(a_q^{m-1})^{\mathrm{T}} \\ b^m(k+1)=b^m(k)-\frac{\alpha}{Q}\sum_{q=1}^Q\delta_q^m \end{matrix}\right.\]
我们这里还是观测权系数$w_{1,1}^1$和$w_{1,1}^2$的批处理学习。下面给出程序：
%初始化已知最优值
W1 = [10; 10];
b1 = [-5;5];
W2 = [1 1];
b2 = [-1];
P = -2:0.1:2;%输入41个值
[R,Q] = size(P);%批处理输入向量
A1 = 1./(1+exp(-(W1*P+b1*ones(1,Q))));
T = 1./(1+exp(-(W2*A1+b2*ones(1,Q))));%真实目标值
lr = 3.5;%学习率
ep = 1000;%学习次数
%*****W1(1,1)，W2(1,1)两个参数是可学习的，其他参数保持不动*****
W1(1,1) = -4;%给定初值
W2(1,1) = -4;
A1 = 1./(1+exp(-(W1*P+b1*ones(1,Q))));
A2 = 1./(1+exp(-(W2*A1+b2*ones(1,Q))));
E = T-A2;%初始的误差向量
SSE0 = sum(sum(E.*E));
xx = [W1(1,1) zeros(1,ep)];
yy = [W2(1,1) zeros(1,ep)];%记录权值的变化
for i=2:(ep+1)
SSE = sum(sum(E.*E));
D2 = A2.*(1-A2).*E;
D1 = A1.*(1-A1).*(W2'*D2);
dW1 = D1*P'*lr;
db1 = D1*ones(Q,1)*lr;
dW2 = D2*A1'*lr;
db2 = D2*ones(Q,1)*lr;
newx = W1(1,1) + dW1(1,1); W1(1,1) = newx; xx(i) = newx;
newy = W2(1,1) + dW2(1,1); W2(1,1) = newy; yy(i) = newy;
A1 = 1./(1+exp(-(W1*P+b1*ones(1,Q))));
A2 = 1./(1+exp(-(W2*A1+b2*ones(1,Q))));
E = T-A2;
ee(i) = sum(sum(E.*E));
end
figure
plot(ee(2:end),'r-o')
set(gca,'FontName','Times New Roman','FontSize',20)
xlabel('epoch times', 'fontsize', 20)%x坐标，设置坐标轴文字大小
ylabel('error square', 'fontsize', 20)%y坐标，设置坐标轴文字大小
figure
W1 =[10;10];
b1 = [-5;5];
W2 = [1 1];
b2 = -1;%正确的权向量和偏置
P=[-2:0.1:2];%输入数据
for k=1:length(P)
p=P(k);
%前向过程
a1 = 1./(1+exp(-(W1*p+b1)));
a2 = 1./(1+exp(-(W2*a1+b2)));
T(k)=a2;
end
%设变量为 W1_{1,1}, W2_{1,1}
[X,Y]=meshgrid(-5:0.5:15);
Z=arrayfun(@(x,y) fcurve(x,y),X,Y);
contour(X,Y,Z,50)
set(gca,'FontName','Times New Roman','FontSize',20)
xlabel('W^1_{1,1}', 'fontsize', 20)%x坐标，设置坐标轴文字大小
ylabel('W^2_{1,1}', 'fontsize', 20)%y坐标，设置坐标轴文字大小
hold on
plot(xx,yy,'r-o')
figure
surf(X,Y,Z)
hold on
plot3(xx(2:end),yy(2:end),ee(2:end),'o-r')
set(gca,'FontName','Times New Roman','FontSize',20)
xlabel('W^1_{1,1}', 'fontsize', 20)%x坐标，设置坐标轴文字大小
ylabel('W^2_{1,1}', 'fontsize', 20)%y坐标，设置坐标轴文字大小
function y=fcurve(w1,w2)
W1 =[10;10];
b1 = [-5;5];
W2 = [1 1];
b2 = -1;
P=[-2:0.1:2];
for k=1:length(P)
p=P(k);
a1 = 1./(1+exp(-(W1*p+b1)));
a2 = 1./(1+exp(-(W2*a1+b2)));
T(k)=a2;
a1 = 1./(1+exp(-([w1;10]*p+b1)));
a2 = 1./(1+exp(-([w2 1]*a1+b2)));
Tw(k)=a2;
end
y=sum((T-Tw).^2);
end
图7.6 对应权系数$w_{1,1}^1$和$w_{1,1}^2$的平方误差学习曲线，误差曲面的等高线和误差曲面中的空间学习曲线
很明显地上面学习曲线在平台区域，总是停滞不前，最后又学习加快，进入坡度平缓的凹槽中。这对于学习率的要求不能太大，否则后面会出现振荡，大家可以改变一下学习率试试。比如学习率变为6，得到如下的曲线
图7.6 对应权系数$w_{1,1}^1$和$w_{1,1}^2$的平方误差学习曲线（大的学习率）
§7.3 冲量改进法
为了提高收敛性，那么可以平滑轨迹，当算法偏离，在凹槽内来回震荡，我们可以采用冲量法。再次回忆权更新公式
\[\left\{\begin{matrix} W^m(k+1)=W^m(k)-\frac{\alpha}{Q}\sum_{q=1}^Q\delta_q^m(a_q^{m-1})^{\mathrm{T}} \\ b^m(k+1)=b^m(k)-\frac{\alpha}{Q}\sum_{q=1}^Q\delta_q^m \end{matrix}\right.\]
\[\left\{\begin{matrix} \Delta W^m(k)=W^m(k+1)-W^m(k)=-\frac{\alpha}{Q}\sum_{q=1}^Q\delta_q^m(a_q^{m-1})^{\mathrm{T}} \\ \Delta b^m(k)=b^m(k+1)-b^m(k)=-\frac{\alpha}{Q}\sum_{q=1}^Q\delta_q^m \end{matrix}\right.\]
引入冲量系数$\gamma$之后，我们改进为
\[\left\{\begin{matrix} \Delta W^m(k)=\gamma \Delta W^m(k)-(1-\gamma)\frac{\alpha}{Q}\sum_{q=1}^Q\delta_q^m(a_q^{m-1})^{\mathrm{T}} \\ \Delta b^m(k)=\gamma \Delta b^m(k)-(1-\gamma)\frac{\alpha}{Q}\sum_{q=1}^Q\delta_q^m \end{matrix}\right.\]
其特点是当轨迹保持在一致的方向上时，会加速收敛。冲量系数$\gamma \geq 0$越大，轨迹保持原有方向的趋势越大。对于上面的问题，我们选用冲量方法，那么得到如下的结果，可以看出振荡消失了，冲量起到了很好的镇定作用。
图7.7 对应权系数$w_{1,1}^1$和$w_{1,1}^2$的平方误差学习曲线（大的学习率+冲量）
%初始化已知最优值
W1 = [10; 10];
b1 = [-5;5];
W2 = [1 1];
b2 = [-1];
P = -2:0.1:2;%输入41个值
[R,Q] = size(P);%批处理输入向量
A1 = 1./(1+exp(-(W1*P+b1*ones(1,Q))));
T = 1./(1+exp(-(W2*A1+b2*ones(1,Q))));%真实目标值
lr = 6;%学习率
mc=0.8;%冲量
ep = 1000;%学习次数
%*****W1(1,1)，W2(1,1)两个参数是可学习的，其他参数保持不动*****
W1(1,1) = -4;%给定初值
W2(1,1) = -4;
A1 = 1./(1+exp(-(W1*P+b1*ones(1,Q))));
A2 = 1./(1+exp(-(W2*A1+b2*ones(1,Q))));
E = T-A2;%初始的误差向量
SSE0 = sum(sum(E.*E));
xx = [W1(1,1) zeros(1,ep)];
yy = [W2(1,1) zeros(1,ep)];%记录权值的变化
dW1 = 0;
db1 = 0;
dW2 = 0;
db2 = 0;
for i=2:(ep+1)
SSE = sum(sum(E.*E));
D2 = A2.*(1-A2).*E;
D1 = A1.*(1-A1).*(W2'*D2);
dW1 = mc*dW1 + (1-mc)*D1*P'*lr;
db1 = mc*db1 + (1-mc)*D1*ones(Q,1)*lr;
dW2 = mc*dW2 + (1-mc)*D2*A1'*lr;
db2 = mc*db2 + (1-mc)*D2*ones(Q,1)*lr;
newx = W1(1,1) + dW1(1,1); W1(1,1) = newx; xx(i) = newx;
newy = W2(1,1) + dW2(1,1); W2(1,1) = newy; yy(i) = newy;
A1 = 1./(1+exp(-(W1*P+b1*ones(1,Q))));
A2 = 1./(1+exp(-(W2*A1+b2*ones(1,Q))));
E = T-A2;
ee(i) = sum(sum(E.*E));
end
figure
plot(ee(2:end),'r-o')
set(gca,'FontName','Times New Roman','FontSize',20)
xlabel('epoch times', 'fontsize', 20)%x坐标，设置坐标轴文字大小
ylabel('error square', 'fontsize', 20)%y坐标，设置坐标轴文字大小
figure
W1 =[10;10];
b1 = [-5;5];
W2 = [1 1];
b2 = -1;%正确的权向量和偏置
%------------
P=[-2:0.1:2];%输入数据
%目标值
for k=1:length(P)
p=P(k);%取一个输入值，可以按顺序取，也可以随机取
%前向过程
a1 = 1./(1+exp(-(W1*p+b1)));
a2 = 1./(1+exp(-(W2*a1+b2)));
T(k)=a2;
end
%设变量为 W1_{1,1}, W2_{1,1}
[X,Y]=meshgrid(-5:0.5:15);
Z=arrayfun(@(x,y) fcurve(x,y),X,Y);
contour(X,Y,Z,50)
set(gca,'FontName','Times New Roman','FontSize',20)
xlabel('W^1_{1,1}', 'fontsize', 20)%x坐标，设置坐标轴文字大小
ylabel('W^2_{1,1}', 'fontsize', 20)%y坐标，设置坐标轴文字大小
hold on
plot(xx,yy,'r-o')
figure
surf(X,Y,Z)
hold on
plot3(xx(2:end),yy(2:end),ee(2:end),'o-r')
set(gca,'FontName','Times New Roman','FontSize',20)
xlabel('W^1_{1,1}', 'fontsize', 20)%x坐标，设置坐标轴文字大小
ylabel('W^2_{1,1}', 'fontsize', 20)%y坐标，设置坐标轴文字大小
function y=fcurve(w1,w2)
W1 =[10;10];
b1 = [-5;5];
W2 = [1 1];
b2 = -1;
P=[-2:0.1:2];
for k=1:length(P)
p=P(k);%取一个输入值，可以按顺序取，也可以随机取
%前向过程
a1 = 1./(1+exp(-(W1*p+b1)));
a2 = 1./(1+exp(-(W2*a1+b2)));
T(k)=a2;
a1 = 1./(1+exp(-([w1;10]*p+b1)));
a2 = 1./(1+exp(-([w2 1]*a1+b2)));
Tw(k)=a2;
end
y=sum((T-Tw).^2);
end
§7.3 数值优化技术-共轭梯度法
回忆二次函数优化的共轭算法步骤
-
给定初始点$x_0$，选择第一个搜索下降方向$p_0=-g_0$
\[g_k=\triangledown F(x)\Bigl|_{x=x_k}\]
-
$a_k=\frac{g_k^{\mathrm{T}}g_k}{g_k^{\mathrm{T}}Ag_k}$, $x_{k+1}=x_k+\alpha_kp_k$
-
$\beta_k=\frac{g_{k+1}^{\mathrm{T}}g_{k+1}}{g_k^{\mathrm{T}}g_k}$
-
$p_{k+1}=-g_{k+1}+\beta_kp_k$
但是因为神经网络的性能函数虽然叫平方误差，但是不是权系数的二次函数，所以不能直接用共轭算法，第二步中要改变为线性搜索。通用过程是特定方向上找到函数的最小值，包括两个步骤：区间定位和区间缩小。区间定位是找到局部极小值的初始区域，区间缩小是缩短初始区域的范围，直到最小值满足精度。
区间定位
区间定位可以用函数比较法来进行，比如上图中$f(x_k+\epsilon p_k)$，$f(x_k+2\epsilon p_k)$，$f(x_k+4\epsilon p_k)$，$f(x_k+8\epsilon p_k)$这样去找点，如果函数连续两次评估都增长，则过程停止，比如上面我们知道在$[a_3,b_3]$ 和$[a_4,b_4]$两个区间内，但是哪一个不知道，因为存在两种情况
上面a-b-c-d直接最小值到底在哪个区间难以用三个函数值进行判定，所以一般采用一些算法不断缩小区间，其中一种为黄金分割搜索。
这里为了可以用一个简单二次函数举例
**************************************
function [i,s,phis,ds,dphi,G]=golds(phi,a,b,epsilon,delta)
%功能: 黄金分割法精确线搜索
%输入: phi是目标函数, a, b是搜索区间的两个端点,
% epsilon, delta分别是自变量和函数值的容许误差
%输出: i是迭代次数, s, phis分别是近似极小点和极小值,
% ds, dphi分别是s和phis的误差限, G是ix4矩阵,
% 其第j行分别是a,p,q,b的第j次迭代值[aj,pj,qj,bj]
t=(sqrt(5)-1)/2; h=b-a;
phia=feval(phi,a); phib=feval(phi,b);
c=a+(1-t)*h; d=a+t*h;
phic=feval(phi,c); phid=feval(phi,d);
i=1; G(i,:)=[a, c, d, b];
while(abs(phib-phia)>delta)|(h>epsilon)
if(phic<=phid)%比较c,d点处函数值，然后bàd, dàc, a保持不动
b=d; phib=phid; d=c; phid=phic;
h=b-a; c=a+(1-t)*h; phic=feval(phi,c); %更新c
else%比较c,d点处函数值，然后aàc, càd, b保持不动
a=c; phia=phic; c=d; phic=phid;
h=b-a; d=a+t*h; phid=feval(phi,d); %更新d
end
i=i+1; G(i,:)=[a, c, d, b];
end
if(phic<=phid)
s=c; phis=phic;
else
s=d; phis=phid;
end
ds=abs(b-a); dphi=abs(phib-phia);
*********************************************
比如我们求函数
在给定的范围为[a=0, b=10]内去找，分别是自变量和函数值的容许误差都是0.1
**********************************************
[i,s,phis,ds,dphi,G]=golds(@(x) (x-1).^2,0,10,1e-1,1e-1)
G =
0 3.8197 6.1803 10.0000
0 2.3607 3.8197 6.1803
0 1.4590 2.3607 3.8197
0 0.9017 1.4590 2.3607
0 0.5573 0.9017 1.4590
0.5573 0.9017 1.1146 1.4590
0.5573 0.7701 0.9017 1.1146
0.7701 0.9017 0.9830 1.1146
0.9017 0.9830 1.0333 1.1146
0.9017 0.9519 0.9830 1.0333
0.9519 0.9830 1.0022 1.0333
a c d b
**********************************************
从运行结果可以看出，最终区间能够缩减到[0.9830,1.0022]。因此数值的共轭梯度方法就是通过选取初始的方向，然后通过区间试探和缩小。然而对于简单的二次函数，n个参数，最多n步可以收敛到最小值，神经网络性能搜索n步之后采用什么方向？有很多提法，最简单的是n步之后重新选为最速下降方向
图7.8 权系数$w_{1,1}^1$和$w_{1,1}^2$的学习曲线
（共轭梯度+线性搜索红色，动量法蓝色）
注意的是虽然上面仅仅5步就到了最低点，但是每一步都需要迭代寻找区间，所以每步都用了很多计算量，但是即便如此，共轭方法还是被证明为多层神经网络的最快训练算法之一。
% INITIAL VALUES
W1 = [10; 10];
b1 = [-5;5];
W2 = [1 1];
b2 = [-1];
P = -2:0.1:2;
[R,Q] = size(P);
A1 = nndlogsig(W1*P+b1*ones(1,Q));%后面定义了就是1/(1+exp(-x))
T = nndlogsig(W2*A1+b2*ones(1,Q));
%% training
W1(1,1) = -4;%给定一对初值
W2(1,1) = -4;
A1 = nndlogsig(W1*P+b1*ones(1,Q));
A2 = nndlogsig(W2*A1+b2*ones(1,Q));
E = T-A2;
fa=sum(sum(E.*E));%初始点平方误差函数
D2 = A2.*(1-A2).*E;%敏感度
D1 = A1.*(1-A1).*(W2'*D2);
gW1 = D1*P';% W1梯度
gb1 = D1*ones(Q,1);
gW2 = D2*A1'; % W2梯度
gb2 = D2*ones(Q,1);
nrmo = gW1(1,1)^2 + gW2(1,1)^2;
% NORM OF GRADIENT
nrmrt=sqrt(nrmo);%梯度
% 初始化方向
dW1old=gW1;db1old=gb1;dW2old=gW2;db2old=gb2;
dW1=gW1/nrmrt;db1=gb1/nrmrt;dW2=gW2/nrmrt;db2=gb2/nrmrt;
% ASSIGN PARAMETERS
tau=0.618;
tau1=1-tau;
scaletol=20;
delta=0.32;
delta1=.03;
tol=delta1/scaletol;
scale=2.0;
bmax=26;%区间[a,b]最大区间值
n=2; %number of steps before reset
% MAIN LOOP
max_epoch = 15;
xx = [W1(1,1) zeros(1,max_epoch)];
yy = [W2(1,1) zeros(1,max_epoch)];
for i=2:(max_epoch+1)
% INITIALIZE A
a=0;
aold=0;
b=delta;
faold=fa;
% CALCULATE INITIAL SSE
W1t = W1; b1t = b1;
W2t = W2; b2t = b2;
newx = W1(1,1) + b*dW1(1,1); W1t(1,1) = newx; % 梯度上走了一步
newy = W2(1,1) + b*dW2(1,1); W2t(1,1) = newy;
fb = sumsqr(T - nndlogsig(W2t*nndlogsig(W1t*P+b1t*ones(1,Q))+b2t*ones(1,Q)));
% FIND INITIAL INTERVAL WHERE SSE MINIMUM OCCURS
while (fa>fb)&(b<bmax)%如果fb试探点函数目标值小了那么继续扩大区间找
aold=a;
faold=fa;
fa=fb;
a=b;
b=scale*b;
newx = W1(1,1) + b*dW1(1,1); W1t(1,1) = newx;
newy = W2(1,1) + b*dW2(1,1); W2t(1,1) = newy;
fb = sumsqr(T - nndlogsig(W2t*nndlogsig(W1t*P+b1t*ones(1,Q))+b2t*ones(1,Q)));
end
a=aold;
fa=faold;%直到试探出来一个区间[a,b]满足fa<fb
% INITIALIZE C AND 进一步缩小区间范围找到局部最小值
c=a+tau1*(b-a);%黄金分割点
newx = W1(1,1) + c*dW1(1,1); W1t(1,1) = newx;
newy = W2(1,1) + c*dW2(1,1); W2t(1,1) = newy;
fc = sumsqr(T - nndlogsig(W2t*nndlogsig(W1t*P+b1t*ones(1,Q))+b2t*ones(1,Q)));
d=b-tau1*(b-a);
newx = W1(1,1) + d*dW1(1,1); W1t(1,1) = newx;
newy = W2(1,1) + d*dW2(1,1); W2t(1,1) = newy;
fd = sumsqr(T - nndlogsig(W2t*nndlogsig(W1t*P+b1t*ones(1,Q))+b2t*ones(1,Q)));
% MINIMIZE ALONG A LINE
k=0;
while (b-a)>tol %一直寻找到区间小于一个允许度
if ( (fc<fd)&(fb>=min([fa fc fd])) ) | fa<min([fb fc fd])
b=d; d=c; fb=fd;
c=a+tau1*(b-a);
fd=fc;
newx = W1(1,1) + c*dW1(1,1); W1t(1,1) = newx;
newy = W2(1,1) + c*dW2(1,1); W2t(1,1) = newy;
fc = sumsqr(T - nndlogsig(W2t*nndlogsig(W1t*P+b1t*ones(1,Q))+b2t*ones(1,Q)));
else
a=c; c=d; fa=fc;
d=b-tau1*(b-a);
fc=fd;
newx = W1(1,1) + d*dW1(1,1); W1t(1,1) = newx;
newy = W2(1,1) + d*dW2(1,1); W2t(1,1) = newy;
fd = sumsqr(T - nndlogsig(W2t*nndlogsig(W1t*P+b1t*ones(1,Q))+b2t*ones(1,Q)));
end
end
% UPDATE VARIABLES
newx = W1(1,1) + a*dW1(1,1); W1(1,1) = newx;
newy = W2(1,1) + a*dW2(1,1); W2(1,1) = newy;
xx(i) = newx;
yy(i) = newy;
% CALCULATE GRADIENT
A1 = nndlogsig(W1*P+b1*ones(1,Q));
A2 = nndlogsig(W2*A1+b2*ones(1,Q));
E = T-A2;
D2 = A2.*(1-A2).*E;
D1 = A1.*(1-A1).*(W2'*D2);
gW1 = D1*P';
gb1 = D1*ones(Q,1);
gW2 = D2*A1';
gb2 = D2*ones(Q,1);
% NORM SQUARE OF GRADIENT
nrmn = gW1(1,1)^2 + gW2(1,1)^2;
% CALCULATE DIRECTION
if rem(i,n)==0 %如果i是n的整数倍则重新回到梯度方向
Z=0;
else
Z=nrmn/nrmo;
end
% CALCULATE NEW DIRECTIONS
dW1new = gW1 + dW1old*Z; db1new = gb1 + db1old*Z;
dW2new = gW2 + dW2old*Z; db2new = gb2 + db2old*Z;
% SAVE NEW DIRECTIONS
dW1old = dW1new; db1old = db1new;
dW2old = dW2new; db2old = db2new;
nrmo=nrmn;
%NORMALIZE DIRECTIONS
nrm = sqrt(dW1new(1,1)^2 + dW2new(1,1)^2);
dW1=dW1new/nrm;db1=db1new/nrm;dW2=dW2new/nrm;db2=db2new/nrm;
end
%figure绘图
W1 =[10;10];
b1 = [-5;5];
W2 = [1 1];
b2 = -1;%正确的权向量和偏置
%------------
P=[-2:0.1:2];%输入数据
%目标值
for k=1:length(P)
p=P(k);%取一个输入值，可以按顺序取，也可以随机取
%前向过程
a1 = 1./(1+exp(-(W1*p+b1)));
a2 = 1./(1+exp(-(W2*a1+b2)));
T(k)=a2;
end
%设变量为 W1_{1,1}, W2_{1,1}
[X,Y]=meshgrid(-5:0.5:15);
Z=arrayfun(@(x,y) fcurve(x,y),X,Y);
contour(X,Y,Z,50)
set(gca,'FontName','Times New Roman','FontSize',20)
xlabel('W^1_{1,1}', 'fontsize', 20)%x坐标，设置坐标轴文字大小
ylabel('W^2_{1,1}', 'fontsize', 20)%y坐标，设置坐标轴文字大小
hold on
plot(xx,yy,'r-o')
function y=fcurve(w1,w2)
W1 =[10;10];
b1 = [-5;5];
W2 = [1 1];
b2 = -1;
P=[-2:0.1:2];
for k=1:length(P)
p=P(k);%取一个输入值，可以按顺序取，也可以随机取
%前向过程
a1 = 1./(1+exp(-(W1*p+b1)));
a2 = 1./(1+exp(-(W2*a1+b2)));
T(k)=a2;
a1 = 1./(1+exp(-([w1;10]*p+b1)));
a2 = 1./(1+exp(-([w2 1]*a1+b2)));
Tw(k)=a2;
end
y=sum((T-Tw).^2);
end
function y=nndlogsig(x)
y=1./(1+exp(-x));
end
§7.4 数值优化技术- Levenberg-Marquardt
LM算法是牛顿方法的变种，非常适合均方误差为性能指标的神经网络训练。我们先回顾一下牛顿法优化函数$f(x)$。利用二阶泰勒展开得出
\[f(x_{k+1})\approx f(x_k)+g_k^{\mathrm{T}}\Delta x_k+\frac{1}{2}\Delta x_k^{\mathrm{T}}A_k\Delta x_k\]
\[g_k=\triangledown f(x)^{\mathrm{T}}\Bigl|_{x=x_k},A_k=\triangledown ^2f(x)^{\mathrm{T}}\Bigl|_{x=x_k}\]
那么该变量如何使得$f(x_{k+1})$最小？对上面展开求导数
\[g_k+A_k\Delta x_k=0\]
我们直接得到
\[\Delta x_k=-A_k^{-1}g_k\]
因此牛顿法的更新规则定义为
\[x_{k+1}=x_k-A_k^{-1}g_k\]
但是这种方法局限性在于因为它不能分辨极小、极大和鞍点。
特别地对于神经网络的均方误差函数
\[F(x)=e^{\mathrm{T}}(x)e(x)=\sum_{i=1}^Ne_i^2(x)\]
这样的多元函数，其梯度的第j个元素
\[[\triangledown F(x)]_j=\frac{\partial F(x)}{\partial x_j}=2\sum_{i=1}^Ne_i(x)\frac{\partial e_i(x)}{\partial x_j}\]
就是Jacobi矩阵形式
\[\triangledown F(x)=2J^{\mathrm{T}}(x)e(x)\]
\[J(x)=\begin{bmatrix} \frac{\partial e_1(x)}{\partial x_1} & \frac{\partial e_1(x)}{\partial x_2} & \cdots & \frac{\partial e_1(x)}{\partial x_N}\\ \frac{\partial e_2(x)}{\partial x_1} & \frac{\partial e_2(x)}{\partial x_2} & \cdots & \frac{\partial e_2(x)}{\partial x_N}\\ \vdots & \vdots & \ddots & \vdots\\ \frac{\partial e_N(x)}{\partial x_1} & \frac{\partial e_N(x)}{\partial x_2} & \cdots & \frac{\partial e_N(x)}{\partial x_N} \end{bmatrix}\]
对于特殊的均方误差函数$F(x)$，其Hessian矩阵的第k,j个元素为
\[[\triangledown ^2F(x)]_{k,j}=\frac{\partial ^2F(x)}{\partial x_k\partial x_j}=2\sum_{i=1}^N\left \{ \frac{\partial e_i(x)}{\partial x_k}\frac{\partial e_i(x)}{\partial x_j}+e_i(x)\frac{\partial ^2e_i(x)}{\partial x_k\partial x_j} \right \}\]
\[\triangledown ^2F(x)=2J^{\mathrm{T}}(x)J(x)+2S(x)\]
\[S(x)=2\sum_{i=1}^{N}\left \{ e_i(x)\frac{\partial ^2e_i(x)}{\partial x_k \partial x_j} \right \}\]
Hessian矩阵计算量太大，现实中经常省去$S(x)$，认为其值很小，因此Hessian矩阵近似为
\[\triangledown ^2F(x)\approx 2J^{\mathrm{T}}(x)J(x)\]
因此更新为
\[x_{k+1}=x_k-A_k^{-1}g_k=x_k-[\triangledown ^2F(x)]^{-1}\triangledown F(x)\]
\[=x_k-[J^{\mathrm{T}}(x)J(x)]^{-1}J^{\mathrm{T}}(x)e(x)\Bigl|_{x=x_k}\]
也称为高斯牛顿法，省去了计算二次导数。但是这种算法经常遇到Hessian矩阵$\triangledown ^2F(x)$不可逆，就是$\triangledown ^2F(x)$可能是半正定的，有特征值接近零。因此Levenberg-Marquardt这样进行了改进，构造一个新的矩阵
\[G=\triangledown ^2F(x)+\mu I,\mu>0\]
这样做的目的很明显，就是把接近零的特征值变成非零。$\triangledown ^2F(x)$的每个特征为$\lambda_i$，$G$的特征值为$\lambda_i+\mu$，这样就变成正定矩阵了。Levenberg-Marquardt算法就是
\[x_{k+1}=x_k-[J^{\mathrm{T}}(x)J(x)+\mu I]^{-1}J^{\mathrm{T}}(x)e(x)\Bigl|_{x=x_k}\]
这个算法复杂度比较大，需要存储逆阵，虽然是最快的BP算法。这里程序仅仅供参考。
图7.8 权系数$w_{1,1}^1$和$w_{1,1}^2$的学习曲线（LM算法）
******
% INITIAL VALUES
W1 = [10; 10];
b1 = [-5;5];
W2 = [1 1];
b2 = [-1];
P = -2:0.1:2;
[R,Q] = size(P);
A1 = nndlogsig(W1*P+b1*ones(1,Q));
T = nndlogsig(W2*A1+b2*ones(1,Q));
%%%
W1(1,1) = -4;
W2(1,1) = -4;
mu_initial =0.1;
v = 2;
A1 = nndlogsig(W1*P+b1*ones(1,Q));
A2 = nndlogsig(W2*A1+b2*ones(1,Q));
E1 = T-A2;
f1 = sum(sum(E1.*E1));
% DEFINE SIZES
[R,Q] = size(P);
[S2,Q] = size(T);
S1 = 2;
RS = S1*R; RS1 = RS+1; RSS = RS + S1; RSS1 = RSS + 1;
RSS2 = RSS + S1*S2; RSS3 = RSS2 + 1; RSS4 = RSS2 + S2;
% ASSIGN PARAMETERS
disp_freq = 1;
max_epoch = 11;
err_goal = 0.0000001;
maxmu=1e10;
mingrad=.000002;
mu=mu_initial;
ii=eye(2);
meu=zeros(max_epoch,1);
mer=meu;grad=meu;
xx = [W1(1,1) zeros(1,max_epoch)];
yy = [W2(1,1) zeros(1,max_epoch)];
for k=2:(max_epoch+1)
% INITIALIZE A
mu=mu/v;
mer(k)=f1;
meu(k)=mu;
tst=1;
% FIND JACOBIAN
A1 = kron(A1,ones(1,S2));
D2 = nnmdlog(A2);
D1 = nnmdlog(A1,D2,W2);
jac1 = nnlmarq(kron(P,ones(1,S2)),D1);
jac2 = nnlmarq(A1,D2);
jac=[jac1,D1',jac2,D2'];
% PULL OUT APPROPRIATE TERMS
jac = [jac(:,1) jac(:,5)];
% CHECK THE MAGNITUDE OF THE GRADIENT
E1=E1(:);
je=jac'*E1;
grad(k)=norm(je);
if grad(k)<mingrad,
mer=mer(1:k);
meu=meu(1:k);
grad=grad(1:k);
disp('Training has stopped.')
disp('Local minimum reached. Gradient is close to zero.')
fprintf('Magnitude of gradient = %g.\n',grad(k));
break
end
% INNER LOOP, INCREASE mu UNTIL THE ERRORS ARE REDUCED
jj=jac'*jac;
while tst>0,
dw=-(jj+ii*mu)\je;
W1n=W1;b1n=b1;W2n=W2;b2n=b2;
% UPDATE VARIABLES
newx = W1(1,1) + dw(1); W1n(1,1) = newx;
newy = W2(1,1) + dw(2); W2n(1,1) = newy;
A1 = nndlogsig(W1n*P+b1n*ones(1,Q));
A2 = nndlogsig(W2n*A1+b2n*ones(1,Q));
E2 = T-A2;
f2=sum(sum(E2.*E2));
if f2>=f1,
mu=mu*v;
% TEST FOR MAXIMUM mu
if (mu > maxmu),
mer=mer(1:k);
meu=[meu(1:k);mu];
grad=grad(1:k);
disp('Maximum mu exceeded.')
fprintf('mu = %g.\n',mu);
fprintf('Maximum allowable mu = %g.\n',maxmu);
break;
end
else
tst=0;
end
end
% TEST IF THE ERROR REACHES THE ERROR GOAL
if f2<=err_goal,
f1=f2;
W1=W1n;b1=b1n;W2=W2n;b2=b2n;
mer=[mer(1:k);f2];
meu=[meu(1:k);mu];
grad=grad(1:k);
disp('Training has stopped. Goal achieved.')
break;
end
if(mu>maxmu),
disp('Maximum mu exceeded.')
fprintf('mu = %g.\n',mu);
fprintf('Maximum allowable mu = %g.\n',maxmu);
break;
end
W1=W1n;b1=b1n;W2=W2n;b2=b2n;E1=E2;f1=f2;
% DISPLAY PROGRESS
if rem(k,disp_freq) == 0
xx(k) = newx;
yy(k) = newy;
end
end
%figure绘图
W1 =[10;10];
b1 = [-5;5];
W2 = [1 1];
b2 = -1;%正确的权向量和偏置
%------------
P=[-2:0.1:2];%输入数据
%目标值
for k=1:length(P)
p=P(k);%取一个输入值，可以按顺序取，也可以随机取
%前向过程
a1 = 1./(1+exp(-(W1*p+b1)));
a2 = 1./(1+exp(-(W2*a1+b2)));
T(k)=a2;
end
%设变量为 W1_{1,1}, W2_{1,1}
[X,Y]=meshgrid(-5:0.5:15);
Z=arrayfun(@(x,y) fcurve(x,y),X,Y);
contour(X,Y,Z,50)
set(gca,'FontName','Times New Roman','FontSize',20)
xlabel('W^1_{1,1}', 'fontsize', 20)%x坐标，设置坐标轴文字大小
ylabel('W^2_{1,1}', 'fontsize', 20)%y坐标，设置坐标轴文字大小
hold on
plot(xx,yy,'r-o')
function y=fcurve(w1,w2)
W1 =[10;10];
b1 = [-5;5];
W2 = [1 1];
b2 = -1;
P=[-2:0.1:2];
for k=1:length(P)
p=P(k);%取一个输入值，可以按顺序取，也可以随机取
%前向过程
a1 = 1./(1+exp(-(W1*p+b1)));
a2 = 1./(1+exp(-(W2*a1+b2)));
T(k)=a2;
a1 = 1./(1+exp(-([w1;10]*p+b1)));
a2 = 1./(1+exp(-([w2 1]*a1+b2)));
Tw(k)=a2;
end
y=sum((T-Tw).^2);
end
function y=nndlogsig(x)
y=1./(1+exp(-x));
end
function d = nnmdlog(a,d,w)
% NNMDLOG Logistic Delta Function for Marquardt.
% Returns the delta values for a layer of
% log-sigmoid neurons for use with Marquardt
% backpropagation.
% (See MDELTALIN,MDELTATAN,LEARN_MARQ,TANSIG)
%
% NNMDLOG(A), A is an SxQ matrix.
% Returns a matrix of delta vectors for an OUTPUT
% layer of log-sigmoid neurons with outputs A
% for the network.
%
% NNMDLOG(A,D,W), A is an S1xQ matrix,
% D is an S2xQ matrix, W is an S2xS1 matrix.
% Returns a matrix of delta vectors for a HIDDEN
% layer of log-sigmoid neurons with outputs A
% which had been passed through a weight matrix W
% to another layer with delta vectors D.
% Copyright 1995-2015 Martin T. Hagan and Howard B. Demuth
% First Version, 8-31-95.
%==================================================================
[s1,q]=size(a);
if nargin == 1
d = -kron(((ones-a).*a),ones(1,s1)).*kron(ones(1,s1),eye(s1));
else
d = ((ones-a).*a).*(w'*d);
end
end
function jac = nnlmarq(p,d)
%NNLMARQ Marquardt Backpropagation Learning Rule
%
% (See PURELIN, LOGSIG, TANSIG)
%
% jac = NNLMARQ(P,D)
% P - RxQ matrix of input vectors.
% D - SxQ matrix of sensitivity vectors.
% Returns:
% jac - a partial jacobian matrix.
% Copyright 1995-2015 Martin T. Hagan and Howard B. Demuth
if nargin ~= 2
error('NNET:nlmarq:Arguments','Wrong number of arguments.');
end
[s,q]=size(d);
[r,q]=size(p);
jac=kron(p',ones(1,s)).*kron(ones(1,r),d');
end
浙公网安备 33010602011771号