本站原创文章，转载请说明来自《老饼讲解-BP神经网络》bp.bbbdata.com

本文是笔者细扒matlab神经网络工具箱traingda的源码后，去除冗余代码，重现的简版traingda代码，代码与traingda的结果完全一致。

通过本代码的学习，可以完全细节的了解用“自适应学习率梯度下降法”求解BP神经网络的实现逻辑。

  01. 代码实现  

  代码简要描述  

本代码主要重现了matlab神经网络工具箱中traingda训练函数训练BP神经网络的逻辑。
并通过一个用例，证明重写后的代码训练出来的权重、阈值与工具箱的训练结果一致。

matlab2014b亲测已跑通：

% 本代码模仿matlab神经网络工具箱的BP神经网络traingda算法,
% 代码主旨用于教学,供大家学习理解matlab工具箱BP神经网络的原理(基于traingda算法)
%  阅读本代码，建议逐行调试
% 转载请说明来自 《老饼讲解神经网络》 bp.bbbdata.com
% ---------特别说明：------------------------
% 本代码仅用于学习，因此不包括工具箱中数据处理部分(即只包括网络构建与训练)，
% 因此，与工具箱对比时，投入的数据需要已经做好归一化处理，并且关闭工具箱的数据分割
function testTrainGda()
% ---------数据生成与参数预设-------------
% 生成输入输出数据
X  = [-1:0.2:1;-1:0.2:1];
y = [sin(X( 1,:)) + X( 2,:);sin(X( 1,:).*X( 1,:)) + 0.5*X( 2,:)];
% 数据归一化
X = 2*(X-repmat(min(X,[],2),1,size(X,2)))./(repmat(max(X,[],2)-min(X,[],2),1,size(X,2)))-1;
y = 2*(y-repmat(min(y,[],2),1,size(y,2)))./(repmat(max(y,[],2)-min(y,[],2),1,size(y,2)))-1;

% 参数预设
hnn     = [8,4];        % 隐节点个数(hideNodeNum)
maxStep = 1000;         % 最大训练步数
goal    = 0.00001;      % 目标误差

%---------调用自写函数进行训练--------------
rand('seed',70);
[W,B,lr] = traingdaBPNet(X,y,hnn,goal,maxStep); % 网络训练
py    = predictBpNet(W,B,X);                   % 网络预测
w12 = W{1,2}                                % 提取网络的权重
w23 = W{2,3}                                % 提取网络的权重
w34 = W{3,4}                                % 提取网络的权重
b2  = B{2};                                 % 提取网络阈值
b3  = B{3};                                 % 提取网络阈值
b4  = B{4};                                 % 提取网络阈值
% -----调用工具箱，与工具箱的结果比较------
rand('seed',70);
net = newff(X,y,hnn,{'tansig','tansig','purelin'},'traingda');


%设置训练参数
net.trainparam.goal        = goal;      % 训练目标
net.trainparam.epochs      = maxStep;   % 最大训练次数.
net.divideParam.trainRatio = 1;         % 全部数据用于训练
net.divideParam.valRatio   = 0;         % 关掉泛化验证数据
net.divideParam.testRatio  = 0;         % 关掉测试数据

% 网络训练
[net,tr,py_tool] = train(net,X,y);   % 训练网络
w12_tool = net.IW{1}                 % 提取网络的权重
w23_tool = net.LW{2,1}               % 提取网络的权重
w34_tool = net.LW{3,2}               % 提取网络的权重

% 与工具箱的差异
maxECompareNet = max([max(abs(py(:)-py_tool(:))),max(abs(w12(:)-w12_tool(:))),max(abs(w23(:)-w23_tool(:))),max(abs(w34(:)-w34_tool(:)))]);
disp(['自写代码与工具箱权重阈值的最大差异:',num2str(maxECompareNet)])
end

function [W,B,lr] = traingdaBPNet(X,y,hnn,goal,maxStep)
%------变量预计算与参数设置-----------
lr        = 0.01;            % 学习率
min_grad  = 1.0e-5;          % 最小梯度
maxE_inc  = 1.04;   %减少学习率的阈值
lr_inc = 1.05;    %学习率升上率
lr_dec = 0.7;     %学习率下降率
%---------初始化WB-------------------
[W,B] = initWB(X,y,hnn);  % 初始化W,B

%---------开始训练--------------------
newW = cell(size( W));  %初始化newW
newB = cell(size( B));  %初始化newB
for t = 1:maxStep
    % 计算当前梯度
    [py,layerVal] = predictBpNet(W,B,X);         % 计算网络的预测值
    [E2,E]        = calMSE(y,py);              % 计算误差
    [gW,gB]       = calGrad(W,layerVal,E);    % 计算梯度
    
    
    %----更新权重阈值-----
    for i = 1:size(W,2)-1
        newW{i,i+1} = W{i,i+1} + lr * gW{i,i+1};%更新梯度
        newB{i+1} = B{i+1} + lr * gB{i+1};%更新阈值
    end
    py = predictBpNet(newW,newB,X);         % 计算网络的预测值
    newE2        = calMSE(y,py);              % 计算误差
    
    if newE2/E2 > maxE_inc   %若果误差上升大于阈值
        lr = lr * lr_dec;    %则降低学习率
        gradLen = 2*min_grad;
    else
        if newE2 < E2        %若果误差减少
            lr = lr * lr_inc;%则增加学习率
        end
        W  = newW;  %更新权值
        B  = newB;  %更新阈值
        gradLen = calGradLen(gW,gB);              % 计算梯度值
    end
    
    %-------检查是否达到退出条件----------
    
    % 如果误差已达要求，或梯度过小,则退出训练
    if E2 < goal   || gradLen <=min_grad
        break;
    end
    
end
end


% ----权值阈值的初始化函数--------------
function [W,B]=initWB(X,y,hnn)
%------------读取网络的基本信息--------
in = size(X,1);    % in: 输入变量的个数(inputNum)
on = size(y,1);    % on: 输出变量个数(outputNum)
nn = [in,hnn,on];  % nn: 各层的节点个数(nodeNum)
ln = size(nn,2);   % ln: 网络层数(layerNum)

% -----计算每层的输入范围和激活范围--------
irMat      = cell(ln,ln);                  % irMat:输入范围矩阵(inputRangeMatrix)，irMat(i,j)代表i层传递给j层的输入范围
irMat{1,2} = [min(X,[],2),max(X,[],2)];    % 第1层传给第2层的输入数据范围
for i = 2:ln-1
    irMat{i,i+1} = repmat([-1,1],nn(i),1);   % 第i层到第i+1层范围为[-1,1]
end
arMat = [inf inf;repmat([-2,2],ln-2,1);inf inf];  % 传递函数的激活范围矩阵(activeRangeMatrix)

%------------逐层初始化权值和阈值------------
W = cell(ln,ln);            % 初始化权重
B = cell(1,ln);            % 初始化阈值
for i = 1 : ln-2
    %----------变量预计算---------------------
    inn   = nn(i);           % inn:输入层神经元个数(inputNodeNum),即第i层的神经元个数
    tnn   = nn(i+1);         % tnn:目标层的神经元个数(targetNodeNum)，即第i+1层的神经元个数
    wMag  = 0.7*tnn^(1/inn); %线性系数
    
    % -------初始化w(按标准输入输出范围)-------------
    w     = 2*rand(tnn,inn)-1;               % 随机初始化权值在[-1,1]之间
    wlen  = sqrt( sum(( w.*w), 2));          % 求出每个目标节点权重长度
    normW = w./repmat( wlen,1,size( w,2));   % w除以权重，即对w单位化
    w     = wMag * normW;                    % w乘以线性系数
    
    % --------初始化b(按标准输入输出范围)----------------
    bIniVal = linspace(-1,1,tnn)';           % 线性生成i+1层的阈值
    bSign   = sign(w(:,1));                  % b的方向要与w中的某一维保持相同
    b       = wMag*bIniVal.*bSign;           % b乘以线性系数
    
    % ---------将wb适配到输入输出范围------------
    %获取本层的输入范围，激活范围
    ir       = irMat{i,i+1};             % 读取输入范围
    irMax    = ir(:,2);                  % 第i层传给i+1层的最大值
    irMin    = ir(:,1);                  % 第i层传给i+1层的最小值
    irL      = (irMax - irMin)/2;        % i层传给i+1层的输入范围长度(向量).
    irCenter = (irMax + irMin)/2;        % i层传给i+1层的输入范围中心(向量).
    
    arMin        = arMat(i+1,1);         % i+1层传递函数的最小激活值
    arMax        = arMat(i+1,2);         % i+1层传递函数的最大激活值
    activeLen    = (arMax - arMin)/2;    % i+1层的激活范围的长度(数值)
    activeCenter = (arMax + arMin)/2;    % i+1层的激活范围的中点(数值)
    
    %适配w,b
    W{i,i+1} = (w ./repmat(irL',tnn,1))*activeLen;                % 权值除以输入范围 再乘以激活范围.
    B{i+1}   = b * activeLen + activeCenter - W{i,i+1} * irCenter;% 阈值扩大到激活范围,并移到激活中心,减去权值的期望输入.
    
end

%最后一个隐层到输出层的传递函数为purelin,不涉及tansig函数，初始值直接随机生成即可
inn        = nn(ln-1);                       % 最后一个隐层的神经元个数
tnn        = nn(ln);                         % 输出层的神经元个数
W{ln-1,ln} = 2 * rand(tnn,inn)-1;            % 随机生成权值范围为[-1,1]
B{ln}      = 2 * rand(tnn,1)-1;              % 随机生成阈值范围为[-1,1]
end


% -------BP的预测函数(返回预测值，和每层的激活值)------
function [y,layerVal] = predictBpNet(W,B,X)%计算各个神经元的值
ln = size(W,2);            % 网络层数
sn = size(X,2);            % 样本个数
layerVal    = cell(1,ln);  % 初始化神经元值cell,

% 计算第1层的值
layerVal{1} = X;           % 第1层的值就是输入

% 逐层计算激活值
for i = 2 : ln-1
    layerVal{i} = tansig(W{i-1,i}*layerVal{i-1}+repmat(B{i},1,sn));
end

%计算最后一层的值
layerVal{ln} = W{ln-1,ln}*layerVal{ln-1}+repmat(B{ln},1,sn);

% 计算网络的输出
y = layerVal{ln} ; % 网络的输出就是最后一层的值
end

% ------误差计算函数-----------------
function [E2,E] = calMSE(y,py)  %计算误差
%用于计算误差，E返回向量，用于计算梯度。E2则是均方误差，用于度量y与y0的误差大小。
e  = y - py;
E  = 2*e/length(y(:));             % 样本的误差梯度，用于后续梯度的计算
E2 = sum(sum(e.*e))/length(y(:));  % 均方差
end


% -------梯度计算函数-------------
function  [gw,gb] = calGrad(W,layerVal,E)
layerNum = size(layerVal,2);    % 网络层数
%-----逐层计算梯度--------
for i = layerNum-1 : -1 : 1   % 计算第i层到i+1层的权值梯度,
    %计算后一层的节点梯度
    if i == layerNum-1        % 最后一层是purelin,需要特殊处理
        gNodeLast = E;
    else
        gNodeLast = W{i+1,i+2}'*gNodeLast.*dtansig(1,layerVal{i+1}); %
    end
    %计算权重阈值梯度
    gw{i,i+1} = gNodeLast* layerVal{i}';
    gb{i+1} = sum(gNodeLast,2);
end
end

% ----梯度长度计算--------------
function gradLen = calGradLen(dwCell,dbCell)
gradList = [];
%将每层的权值和阈值合并到向量中
for i = 1 : size(dwCell)-1
    curdW = dwCell{i,i+1}; %
    gradList =  [gradList;curdW(:);dbCell{i+1}];
end
%计算梯度长度
gradLen = sqrt(sum(sum(gradList.^2)));
end

  02 .traingda与traingd的区别  

traingda与traingd基本是一致的，它只是在traingd的基础上加入了自适应学习率。

自适应学习率的机制主要有两点：
(1) 新解的误差允许上升，但不能上升太多。如果误差上升较多，则降低学习率，用更小的学习率去迭代。
(2) 如果误差减少，则加大学习率，即增大调整步长，通过这样的机制使收敛速度更加快。                       

  03. 代码运行结果解说    

运行结果共三部分

1. 自写代码求得的网络权重与阈值

2. 调用工具箱求得的网络权重与阈值

3. 自写代码与工具箱的结果对比

从运行结果可以看到，自写代码与工具箱的结果一样，说明扒出的逻辑与工具箱的基本一致。

与工具箱结果差异的来源分析

和工具箱还是存在一些不可忽视的误差，
经笔者仔细跟进与分析后，确定误差主要来源于数值问题，
虽然刚开始误差只有10^(-13)，但由于学习率存在阈值判断，
当新旧解的误差比刚好在阈值附近时，自写代码与工具箱改变学习率的方向就会不一样，导致后面误差越来越大。

  04. 代码结构说明  

代码主要包含了4个主要函数和3个辅助函数：

testTrainGda：测试用例主函数，直接运行时就是执行该函数。

1、数据：生成一个2输入2输出的训练数据，
2、用自写的函数构建一个2隐层BP神经网络（使用梯度下降法训练网络）
3、使用工具箱newff训练一个2隐层BP神经网络(使用梯度下降法-traingd训练网络)。
4、比较自写函数与工具箱训练结果是否一致(比较测试数据的输出结果是否一致)

traingdaBPNet：训练主函数。

使用梯度下降法训练一个BP神经网络。

initWB：BP权重阈值初始化函数。

使用nguyen_Widrow法初始化BP的权重和阈值。

predictBpNet：预测主函数。

传入需要预测的X，与网络的权重矩阵，即可得到预测结果。

其它辅助函数：

calMSE：误差计算函数
calGrad：梯度计算函数
gradLen：梯度长度计算