标准BP算法matlab实现,简单易懂

机器学习的实验课要求自编写一份简易的标准BP（神经网络）算法,我用matlab基本实现了，现将自己的思想记录下来，方便自己以后重温.话不多说，让我们进入正题

调用matlab的神经网络算法解决具体问题<参照我另一篇博客>

这里是自己编写的python实现标准BP算法

文章目录

简要介绍
标准BP算法
累积BP算法

简要介绍

神经网络中最基本的成分是神经元模型,先简单提一个最简单的神经元模型M-P神经元模型

该神经元模型有两层：输入层,输出层。
其中输入层包含 n 个神经元,输出层包含 1 个神经元。
输入信号通过连接的权重wiw_iwi到达输出神经元,总接收的输入值与阀值θ\thetaθ进行比较,然后通过"激活函数"处理以产生神经元的输出。这里介绍一个最具典型的Sigmoid函数

Sigmod函数能把较大范围内变化的输入值积压到(0,1)的输出值范围内,有时也被称为"挤压函数"

通常情况下,为了处理更加复杂的问题,在输入层与输出层之间会增添多个隐层.常见的神经网络形如下图的层级结构,每层神经元与下层神经元全互连,同层神经元不进行连接,也不存在跨层连接。其中需要注意的是:

输入层神经元接收外界输入,隐层与输出层神经元对信号进行加工,最终结果由输出层神经元输出.
连接的边都具有相应的权值，隐层和输出层的神经单元上还有阀值和激活函数.

前期介绍完毕，接下来进行误差逆传播(error BackPropagation 简称 BP )算法的介绍

标准BP算法

上图是一个三层神经元模型,其中输入层神经元数目为ddd,隐层神经元数目为qqq,输出层神经元数目为lll,
标准BP算法主要针对一个训练样例（假定输入是Rd\mathbb{R^d}Rd维,输出是Rl\mathbb{R^l}Rl维）更新连接权值和阀值.
这里先对一些符号进行说明：

xix_ixi：传输信号,第i个输入神经元的传入信息
vihv_{ih}vih：第i个输入神经元与第h个隐层神经元的连接权值
αh\alpha_hαh：第h个隐层神经元从输入神经层接收到的信息值,αh=∑i=1dvihxi\alpha_h=\sum_{i=1}^dv_{ih}x_iαh=∑i=1dvihxi
γh\gamma_hγh：第h个隐层神经元的阀值(虽未在图中标识,但实际存在)
bhb_hbh：第h个隐层神经元的输出值,有公式:bh=f1(αh−γh)b_h=f_1(\alpha_h-\gamma_h)bh=f1(αh−γh),f1f_1f1为隐层的激活函数
whjw_{hj}whj：第h个隐层神经元与第j个输出神经元的连接权值
βj\beta_jβj：第j个输出神经元从隐层神经层接收到的信息值,βj=∑h=1qwhjbh\beta_j=\sum_{h=1}^qw_{hj}b_hβj=∑h=1qwhjbh
θj\theta_jθj：第j个输出神经元的阀值(虽未在图中标识,但实际存在)
yjy_jyj：第j个输出层神经元的输出值,有公式:yj=f2(βj−θj)y_j=f_2(\beta_j-\theta_j)yj=f2(βj−θj),f2f_2f2为输出层的激活函数

f1,f2f_1,f_2f1,f2一般取sigmod函数，其次分析上面的参数,我们已知xi,yjx_i,y_jxi,yj,我们需要去求vih,γh,whj,θjv_{ih},\gamma_h,w_{hj},\theta_jvih,γh,whj,θj，

需要求解的参数个数为：(d+l+1)q+l,其中:
输入层到隐层连边权值:d*q；隐层的阀值:q；隐层到输出层的权值：q*l；输出层的阀值：l
总计：d*q+q+q*l+l=(d+l+1)q+l

我们需要去求解一个比较好的连边权值和阀值。初始随机设定连边权值和阀值,通过计算得到的误差反向修正连边权值和阀值。那么误差怎么计算呢,通过实际的输出和BP神经网络得到的yjy_jyj进行比较得到误差.

标准BP算法针对的是一个训练样例(注意是一个样本！！)，假设输入为X={x1,x2,...,xd}∈RdX=\{x_1,x_2,...,x_d\}\in\mathbb{R^d}X={x1,x2,...,xd}∈Rd,输出Y={y1,y2,..,yl}∈RlY=\{y_1,y_2,..,y_l\}\in\mathbb{R^l}Y={y1,y2,..,yl}∈Rl。

通过标准BP算法得到的输出为Y^\hat{Y}Y^={y1^,y2^,...,yl^\hat{y_1},\hat{y_2},...,\hat{y_l}y1^,y2^,...,yl^}那么此时该样本的均方误差为:
E=12∑j=1l(yj^−yj)2E=\frac{1}{2}\sum_{j=1}^l(\hat{y_j}-y_j)^2 E=21j=1∑l(yj^−yj)2
此处前面乘以1/2是人为这样写，方便后面计算,我们通过均方误差对连边权值和阀值进行反向修正

即有：

v←v+Δvv\leftarrow v+\Delta vv←v+Δv
γ←γ+Δγ\gamma\leftarrow \gamma+\Delta \gammaγ←γ+Δγ
w←w+Δww\leftarrow w+\Delta ww←w+Δw
θ←θ+Δθ\theta\leftarrow \theta+\Delta \thetaθ←θ+Δθ

这里主要运用了梯度下降策略,以目标的负梯度方向对参数进行调整,对于误差E，给定学习率η\etaη(一般取0.1)，以隐层到输出层的连边权值www为例子:
Δwhj=−η∂E∂whj\Delta w_{hj}=-\eta \frac{\partial E}{\partial w_{hj}} Δwhj=−η∂whj∂E
我们能发现如下的影响链条whj→βj→yj^→Ew_{hj}\rightarrow \beta_j \rightarrow \hat{y_j} \rightarrow Ewhj→βj→yj^→E
∂E∂whj=∂E∂yj^∗∂yj^∂βj∗∂βj∂whj\frac{\partial E}{\partial w_{hj}}=\frac{\partial E}{\partial \hat{y_j}}*\frac{\partial \hat{y_j}}{\partial \beta_j}*\frac{\partial \beta_j}{\partial w_{hj}} ∂whj∂E=∂yj^∂E∗∂βj∂yj^∗∂whj∂βj

∂E∂yj^=yj^−yj\frac{\partial E}{\partial \hat{y_j}}=\hat{y_j}-y_j∂yj^∂E=yj^−yj
∂yj^∂βj=f2′(βj−θj)=yj^(1−yj^)\frac{\partial \hat{y_j}}{\partial \beta_j}=f'_2(\beta_j-\theta_j)=\hat{y_j}(1-\hat{y_j})∂βj∂yj^=f2′(βj−θj)=yj^(1−yj^)
∂βj∂whj=bh\frac{\partial \beta_j}{\partial w_{hj}}=b_h∂whj∂βj=bh

解释以下上面第二个公式,f2′(βj−θj)=yj^(1−yj^)f'_2(\beta_j-\theta_j)=\hat{y_j}(1-\hat{y_j})f2′(βj−θj)=yj^(1−yj^) 注意：yj=f2(βj−θj)y_j=f_2(\beta_j-\theta_j)yj=f2(βj−θj)
Sigmod函数的性质:
f′(x)=(11+e−x)′=e−x(1+e−x)2=11+e−x∗e−x1+e−x=f(x)(1−f(x))f'(x)=(\frac{1}{1+e^{-x}})'=\frac{e^{-x}}{(1+e^{-x})^2}=\frac{1}{1+e^{-x}}*\frac{e^{-x}}{1+e^{-x}}=f(x)(1-f(x)) f′(x)=(1+e−x1)′=(1+e−x)2e−x=1+e−x1∗1+e−xe−x=f(x)(1−f(x))
所以有:
Δwhj=−η∂E∂whj=−η∗(yj^−yj)∗yj^(1−yj^)∗bh\Delta w_{hj}=-\eta \frac{\partial E}{\partial w_{hj}}=-\eta*(\hat{y_j}-y_j)*\hat{y_j}(1-\hat{y_j})*b_h Δwhj=−η∂whj∂E=−η∗(yj^−yj)∗yj^(1−yj^)∗bh
若令gj=(yj−yj^)∗yj^(1−yj^)g_j=(y_j-\hat{y_j})*\hat{y_j}(1-\hat{y_j})gj=(yj−yj^)∗yj^(1−yj^),上式可化简为:
Δwhj=η∗gj∗bh\Delta w_{hj}=\eta*g_j*b_h Δwhj=η∗gj∗bh

类似可通过梯度下降得到其余三个,即

Δθj=−η∂E∂θj=−ηgj\Delta \theta_j=-\eta \frac{\partial E}{\partial \theta_j}=-\eta g_jΔθj=−η∂θj∂E=−ηgj
Δvih=∂E∂vih=ηehxi\Delta v_{ih}=\frac{\partial E}{\partial v_{ih}}=\eta e_hx_iΔvih=∂vih∂E=ηehxi
Δγh=∂E∂γh−ηeh\Delta \gamma_h=\frac{\partial E}{\partial \gamma_h}-\eta e_hΔγh=∂γh∂E−ηeh

其中
eh=−∂E∂bh∂bh∂αh=−∑j=1l∂E∂βj∂βj∂bhf′(αh−γh)=∑j=1lwhjgjf′(αh−γh)=bh(1−bh)∑j=1lwhjgje_h=-\frac{\partial E}{\partial b_h}\frac{\partial b_h}{\partial \alpha_h}=-\sum_{j=1}^l\frac{\partial E}{\partial \beta_j}\frac{\partial \beta_j}{\partial b_h}f'(\alpha_h-\gamma_h)=\sum_{j=1}^lw_{hj}g_jf'(\alpha_h-\gamma_h)=b_h(1-b_h)\sum_{j=1}^lw_{hj}g_j eh=−∂bh∂E∂αh∂bh=−j=1∑l∂βj∂E∂bh∂βjf′(αh−γh)=j=1∑lwhjgjf′(αh−γh)=bh(1−bh)j=1∑lwhjgj

ok，到这里有关标准BP算法的介绍已经讲完了,附上实现的matlab代码

clear;clc
%标准BP算法
x=[1.24 1.27];y=[1 0];  %初始数据
[v,r,w,h,y_hat]=standard_BP(x,y,3,1e-6);function [v,r,w,h,y_hat]=standard_BP(x,y,q,eps)  %q为隐层单元数目,eps均方误差限q=3;    %隐层单元数目
L=length(y);    %输出单元数目
n=length(x);  %获取数据的维度
v=rand(n,q);  %初始化输入层到隐层的权值
r=rand(1,q);    %初始化隐层的阀值
w=rand(q,L);  %初始化隐层到输出层的权值
h=rand(1,L);    %初始化输出层的阀值
k=0.1;        %学习率
E=1;
while E>eps  A=x*v;                %输入层->隐层,各个隐层单元具有的权值
b=fc_sigmod(A-r);     %经过隐层的激活函数的输出B=b*w;                %隐层->输出层,各个输出层单元具有的权值
y_hat=fc_sigmod(B-h); %经过输出层的激活函数的输出E=0.5*sum((y_hat-y).^2);  %求均方误差%以下对各个系数进行调整
g=y_hat.*(1-y_hat).*(y-y_hat);
e=b.*(1-b).*(w*g')';for i=1:nfor j=1:qv(i,j)=v(i,j)+k*e(j)*x(i);   %输入层->隐层的权值更新end
end
r=r-k*e;         %隐层的阀值更新for i=1:qfor j=1:Lw(i,j)=w(i,j)+k*g(j)*b(i);   %隐层->输出层的权值更新end
end
h=h-k*g;enddisp('输入层到隐层的权值为');v
disp('隐层的阀值为');r
disp('隐层到输出层的权值为');w
disp('输出层的阀值为');hfunction y=fc_sigmod(x)
y=1./(1+exp(-x));
end
end

我输入数据X=[1.24 1.27],输出数据是Y=[1 0],

拟合基本达到了要求,这里的程序不限定样本的输入和输出的维度,可自行调控。

数据集通常是包含多个样本.在标准BP算法中，对每个样本都会进行更新连边权值和阀值,对一组数据进行训练.这里给出更为一般的标准BP算法.

给出数据集，15个样本,输入与输出的维度都为2.

x=[1.24,1.27;1.36,1.74;1.38,1.64;1.38,1.82;1.38,1.90;1.40,1.70;1.48,1.82;1.54,1.82;1.56,2.08;1.14,1.82;1.18,1.96;1.20,1.86;1.26,2.00 ;1.28,2.00;1.30,1.96];
y=[1 0;1 0;1 0;1 0;1 0;1 0;1 0;1 0;1 0;0 1;0 1;0 1;0 1;0 1;0 1];

这里跳出调整参数，我改成了迭代次数,

clear;clc
%标准BP算法
x=[1.24,1.27;1.36,1.74;1.38,1.64;1.38,1.82;1.38,1.90; 1.40,1.70;1.48,1.82;1.54,1.82;1.56,2.08;1.14,1.82;1.18,1.96;1.20,1.86;1.26,2.00 1.28,2.00;1.30,1.96];
y=[1 0;1 0;1 0;1 0;1 0;1 0;1 0;1 0;1 0;0 1;0 1;0 1;0 1;0 1;0 1];  %初始数据[v,r,w,h,y_hat]=standard_BP(x,y,3,1e6);x_k=[1.24 1.80;1.28 1.84;1.40 2.04];function [v,r,w,h,y_hat]=standard_BP(x0,y0,q,N)  %q为隐层单元数目,N,迭代次数%初始化连边权值和阀值
L=size(y0,2);    %输出单元数目
[m,n]=size(x0);  %获取数据的维度
v=rand(n,q);  %初始化输入层到隐层的权值
r=rand(1,q);    %初始化隐层的阀值
w=rand(q,L);  %初始化隐层到输出层的权值
h=rand(1,L);    %初始化输出层的阀值
k=0.05;        %学习率
E=1;iter=1;
a=1;while iter<Nx=x0(a,:);y=y0(a,:);A=x*v;                %输入层->隐层,各个隐层单元具有的权值b=fc_sigmod(A-r);     %经过隐层的激活函数的输出B=b*w;                %隐层->输出层,各个输出层单元具有的权值y_hat(a,:)=fc_sigmod(B-h); %经过输出层的激活函数的输出E=0.5*sum((y_hat(a,:)-y).^2);  %求均方误差%以下对各个系数进行调整g=y_hat(a,:).*(1-y_hat(a,:)).*(y-y_hat(a,:));e=b.*(1-b).*(w*g')';for i=1:nfor j=1:qv(i,j)=v(i,j)+k*e(j)*x(i);   %输入层->隐层的权值更新endendr=r-k*e;         %隐层的阀值更新for i=1:qfor j=1:Lw(i,j)=w(i,j)+k*g(j)*b(i);   %隐层->输出层的权值更新endendh=h-k*g;if a>=ma=a-m;enda=a+1;iter=iter+1;
enddisp('输入层到隐层的权值为');v
disp('隐层的阀值为');r
disp('隐层到输出层的权值为');w
disp('输出层的阀值为');hfunction y=fc_sigmod(x)
y=1./(1+exp(-x));
end
end

调整后，输出的y_hat值达到较高的精度

累积BP算法

因为标准BP算法每次更新都只针对单个样例,参数更新得非常频繁,而且对不同样例进行更新的效果可能出现“抵消”现象,因此，为了达到同样的累积误差极小点,标准BP算法往往需要进行更多次数的迭代.累积BP算法直接针对累积误差最小化,它在读取整个训练集D一遍后才对参数进行更新,其参数更新的频率低得多.
这里给出累积误差公式:
E^=1m∑k=1mEk\hat{E}=\frac{1}{m}\sum_{k=1}^mE_k E^=m1k=1∑mEk
我们用累积误差E^\hat{E}E^去对参数进行调整,举一个例子,比如：

Δwhj=−η∂E^∂whj=−η∑k=1m∂E^∂Ek∂Ek∂whj=−η1m∑k=1m∂Ek∂whj\Delta w_{hj}=-\eta\frac{\partial \hat{E}}{\partial w_{hj}}=-\eta\sum_{k=1}^m\frac{\partial \hat{E}}{\partial E_k}\frac{\partial E_k}{\partial w_{hj}}=-\eta\frac{1}{m}\sum_{k=1}^m\frac{\partial E_k}{\partial w_{hj}}Δwhj=−η∂whj∂E^=−η∑k=1m∂Ek∂E^∂whj∂Ek=−ηm1∑k=1m∂whj∂Ek
计算出m个样本的∂Ek∂whj\frac{\partial E_k}{\partial w_{hj}}∂whj∂Ek后取均值,即可

当然,我们可以对上式进行进一步化简,由标准BP算法，我们有：

Δwhj=ηgjbh\Delta w_{hj}=\eta g_jb_hΔwhj=ηgjbh
Δθj=−η∂E∂θj=−ηgj\Delta \theta_j=-\eta \frac{\partial E}{\partial \theta_j}=-\eta g_jΔθj=−η∂θj∂E=−ηgj
Δvih=ηehxi\Delta v_{ih}=\eta e_hx_iΔvih=ηehxi
Δγh=−ηeh\Delta \gamma_h=-\eta e_hΔγh=−ηeh

我们只需要把对应的g_j,b_h,e_h以及x_i对应位置上通过m个样本计算完后,在其对应位置上取均值即可。

这里附上累积BP算法:

clear;clc
%累积BP算法
x=[1.24,1.27;1.36,1.74;1.38,1.64;1.38,1.82;1.38,1.90; 1.40,1.70;1.48,1.82;1.54,1.82;1.56,2.08;1.14,1.82;1.18,1.96;1.20,1.86;1.26,2.00 1.28,2.00;1.30,1.96];
y=[1 0;1 0;1 0;1 0;1 0;1 0;1 0;1 0;1 0;0 1;0 1;0 1;0 1;0 1;0 1];  %初始数据[v,r,w,h,y_hat]=acc_BP(x,y,3,1e3);function [v,r,w,h,y_hat]=acc_BP(x0,y0,q,N)  %q为隐层单元数目,N为迭代次数%初始化连边权值和阀值
L=size(y0,2);    %输出单元数目
[m,n]=size(x0);  %获取样本个数以及数据维度
v=rand(n,q);  %初始化输入层到隐层的权值
r=rand(1,q);    %初始化隐层的阀值
w=rand(q,L);  %初始化隐层到输出层的权值
h=rand(1,L);    %初始化输出层的阀值
k=0.1;        %学习率iter=1;
while iter<Nfor i=1:mx=x0(i,:);y=y0(i,:);A=x*v;                %输入层->隐层,各个隐层单元具有的权值b(i,:)=fc_sigmod(A-r);     %经过隐层的激活函数的输出B=b(i,:)*w;                %隐层->输出层,各个输出层单元具有的权值
y_hat(i,:)=fc_sigmod(B-h); %经过输出层的激活函数的输出%以下对各个系数进行调整
g(i,:)=y_hat(i,:).*(1-y_hat(i,:)).*(y-y_hat(i,:));
e(i,:)=b(i,:).*(1-b(i,:)).*(w*g(i,:)')';
end%对上述梯度下降策略后的调整至取均值
b_bar=mean(b);
g_bar=mean(g);
e_bar=mean(e);    E=mean(sum((y_hat-y0).^2,2));  %求均方误差for i=1:nfor j=1:qv(i,j)=v(i,j)+k*e_bar(j)*mean(x0(:,i));   %输入层->隐层的权值更新endendr=r-k*e_bar;         %隐层的阀值更新for i=1:qfor j=1:Lw(i,j)=w(i,j)+k*g_bar(j)*b_bar(i);   %隐层->输出层的权值更新endendh=h-k*g_bar;        %输出层的阀值更新iter=iter+1;
end
disp('输入层到隐层的权值为');v
disp('隐层的阀值为');r
disp('隐层到输出层的权值为');w
disp('输出层的阀值为');h
disp('此时误差E为');Efunction y=fc_sigmod(x)
y=1./(1+exp(-x));
end
end

这里调整后求出的y^\hat{y}y^效果不是很好，后期数据收敛了。

我们通过神经网络的内在结构求出各个连边权值和阀值,由有标签的数据，去预测无标签的数据类型。最后，我附上自己编写的预测代码

function y_hat=standard_BP_predict(data,v,r,w,h)
[m,n]=size(data);
for i=1:mx=data(i,:);A=x*v;                %输入层->隐层,各个隐层单元具有的权值b=fc_sigmod(A-r);     %经过隐层的激活函数的输出B=b*w;                %隐层->输出层,各个输出层单元具有的权值y_hat(i,:)=fc_sigmod(B-h); %经过输出层的激活函数的输出
end
function y=fc_sigmod(x)
y=1./(1+exp(-x));
end
end

data是我们需要去预测的数据,得到其类型,而v,r,w,h是前期可通过标准BP算法或者累积BP算法算得。

当然，想要获得精度更准确的数据，还得参考其他大量有关BP算法的文章。本文只是简易的对标准BP算法以及累积BP算法做个实现.