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

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

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

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

简要介绍

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

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

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

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

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

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

标准BP算法

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

• $x_i$：传输信号,第i个输入神经元的传入信息
• $v_{ih}$：第i个输入神经元与第h个隐层神经元的连接权值
• ${\alpha }_h$：第h个隐层神经元从输入神经层接收到的信息值,${\alpha }_h={\sum }_{i=1}^d{v}_{ih}{x}_i$
• ${\gamma }_h$：第h个隐层神经元的阀值(虽未在图中标识,但实际存在)
• $b_h$：第h个隐层神经元的输出值,有公式:$b_h=f_1({\alpha }_h-{\gamma }_h)$,$f_1$为隐层的激活函数
• $w_{hj}$：第h个隐层神经元与第j个输出神经元的连接权值
• ${\beta }_j$：第j个输出神经元从隐层神经层接收到的信息值,${\beta }_j={\sum }_{h=1}^q{w}_{hj}{b}_h$
• ${\theta }_j$：第j个输出神经元的阀值(虽未在图中标识,但实际存在)
• $y_j$：第j个输出层神经元的输出值,有公式:$y_j=f_2({\beta }_j-{\theta }_j)$,$f_2$为输出层的激活函数

$f_1,f_2$一般取sigmod函数，其次分析上面的参数,我们已知$x_i,y_j$,我们需要去求$v_{ih},{\gamma }_h,w_{hj},{\theta }_j$，

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

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

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

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

即有：

• $v\leftarrow v+\Delta v$
• $\gamma \leftarrow \gamma +\Delta \gamma$
• $w\leftarrow w+\Delta w$
• $\theta \leftarrow \theta +\Delta \theta$

这里主要运用了梯度下降策略,以目标的负梯度方向对参数进行调整,对于误差E，给定学习率$\eta$(一般取0.1)，以隐层到输出层的连边权值$w$为例子:
$$\Delta w_{hj}=-\eta \frac{\partial E}{\partial w_{hj}}$$
我们能发现如下的影响链条$w_{hj}\to {\beta }_j\to \widehat{y_j}\to E$
$$\frac{\partial E}{\partial w_{hj}}=\frac{\partial E}{\partial \widehat{y_j}}\ast \frac{\partial \widehat{y_j}}{\partial {\beta }_j}\ast \frac{\partial {\beta }_j}{\partial w_{hj}}$$

• $\frac{\partial E}{\partial \widehat{y_j}}=\widehat{y_j}-y_j$
• $\frac{\partial \widehat{y_j}}{\partial {\beta }_j}=f_2^{\prime }({\beta }_j-{\theta }_j)=\widehat{y_j}(1-\widehat{y_j})$
• $\frac{\partial {\beta }_j}{\partial w_{hj}}=b_h$

解释以下上面第二个公式,$f_2^{\prime }({\beta }_j-{\theta }_j)=\widehat{y_j}(1-\widehat{y_j})$ 注意：$y_j=f_2({\beta }_j-{\theta }_j)$
Sigmod函数的性质:
$$f^{\prime }(x)=(\frac{1}{1+e^{-x}}{)}^{\prime }=\frac{e^{-x}}{(1+e^{-x}{)}^2}=\frac{1}{1+e^{-x}}\ast \frac{e^{-x}}{1+e^{-x}}=f(x)(1-f(x))$$
所以有:
$$\Delta w_{hj}=-\eta \frac{\partial E}{\partial w_{hj}}=-\eta \ast (\widehat{y_j}-y_j)\ast \widehat{y_j}(1-\widehat{y_j})\ast b_h$$
若令$g_j=(y_j-\widehat{y_j})\ast \widehat{y_j}(1-\widehat{y_j})$,上式可化简为:
$$\Delta w_{hj}=\eta \ast g_j\ast b_h$$

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

• $\Delta {\theta }_j=-\eta \frac{\partial E}{\partial {\theta }_j}=-\eta g_j$
• $\Delta v_{ih}=\frac{\partial E}{\partial v_{ih}}=\eta e_h{x}_i$
• $\Delta {\gamma }_h=\frac{\partial E}{\partial {\gamma }_h}-\eta e_h$

其中
$$e_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}^{\prime }({\alpha }_h-{\gamma }_h)=\sum _{j=1}^l{w}_{hj}{g}_j{f}^{\prime }({\alpha }_h-{\gamma }_h)=b_h(1-b_h)\sum _{j=1}^l{w}_{hj}{g}_j$$

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:n
    for j=1:q
        v(i,j)=v(i,j)+k*e(j)*x(i);   %输入层->隐层的权值更新
    end
end
r=r-k*e;         %隐层的阀值更新

for i=1:q
    for j=1:L
        w(i,j)=w(i,j)+k*g(j)*b(i);   %隐层->输出层的权值更新
    end
end
h=h-k*g;

end

disp('输入层到隐层的权值为');v
disp('隐层的阀值为');r
disp('隐层到输出层的权值为');w
disp('输出层的阀值为');h

function 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<N
    x=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:n
        for j=1:q
            v(i,j)=v(i,j)+k*e(j)*x(i);   %输入层->隐层的权值更新
        end
    end
    r=r-k*e;         %隐层的阀值更新
    
    for i=1:q
        for j=1:L
            w(i,j)=w(i,j)+k*g(j)*b(i);   %隐层->输出层的权值更新
        end
    end
    h=h-k*g;
    
    if a>=m
        a=a-m;
    end
    a=a+1;
    iter=iter+1;
end

disp('输入层到隐层的权值为');v
disp('隐层的阀值为');r
disp('隐层到输出层的权值为');w
disp('输出层的阀值为');h

function y=fc_sigmod(x)
y=1./(1+exp(-x));
end
end

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

---

累积BP算法

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

• $\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}}$
  计算出m个样本的$\frac{\partial E_k}{\partial w_{hj}}$后取均值,即可

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

• $\Delta w_{hj}=\eta g_j{b}_h$
• $\Delta {\theta }_j=-\eta \frac{\partial E}{\partial {\theta }_j}=-\eta g_j$
• $\Delta v_{ih}=\eta e_h{x}_i$
• $\Delta {\gamma }_h=-\eta e_h$

我们只需要把对应的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<N

for i=1:m
    x=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:n
        for j=1:q
            v(i,j)=v(i,j)+k*e_bar(j)*mean(x0(:,i));   %输入层->隐层的权值更新
        end
    end
    r=r-k*e_bar;         %隐层的阀值更新
    
    for i=1:q
        for j=1:L
            w(i,j)=w(i,j)+k*g_bar(j)*b_bar(i);   %隐层->输出层的权值更新
        end
    end
    h=h-k*g_bar;        %输出层的阀值更新
    
iter=iter+1;    
end
disp('输入层到隐层的权值为');v
disp('隐层的阀值为');r
disp('隐层到输出层的权值为');w
disp('输出层的阀值为');h
disp('此时误差E为');E


function y=fc_sigmod(x)
y=1./(1+exp(-x));
end
end

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

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

function y_hat=standard_BP_predict(data,v,r,w,h)
[m,n]=size(data);
for i=1:m
    x=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算法做个实现.