“手撕”BP算法——使用MATLAB搭建简单的神经网络（附代码）
之前一直都是直接使用深度学习的框架，但对里面所涉及到的基本算法却没有深入研究。看了吴恩达的机器学习视频之后，决定使用MATLAB实现一个简单的神经网络，深刻体会到只有用代码从头实现一个算法，才会对这个算法理解得更加深刻，也才能真正掌握该算法。 机器学习定义如下：一个程序被认为能从经验E中学习，解决任务T，达到性能P，当且仅当，有了经验E之后，经过度量P的评判，程序在处理T的性能有所提升。 神经网络是机器学习中的一类算法，在训练过程中，神经网络内部不断地调整其内部参数的大小，使得神经网络的输出不断地向标签靠拢。其中，“调整参数大小”这个过程一般是“梯度下降”，在一个多层的神经网络中，要达到“梯度下降”这个目的，就不得不提一下反向传播算法(Backpropagation Algorithm, BP)。 BP算法的核心就是链式求导法则，这里对BP算法的理论推导如下，主要参考[1]-[3]。 如图所示，考虑一个有三层、每层有两个神经元的神经网络： 上图为决策面可视化结果，蓝色为判为(0,1)与(1,0)的区域，红色为判为(0,0)，(1,1)的区域。 上图为修改损失函数为交叉熵时的结果，学习率同样为0.5，可以看出在400次左右便达到收敛，因此损失函数的设计对整个神经网络的参数调整至关重要，至于为什么将损失函数调整为交叉熵后，收敛速度会增加，可自行查阅资料，这里给出一个结论：当输出层采用线性激活函数时，损失函数采用MSE会收敛更快；当输出层采用sigmoid函数时，一般采用交叉熵会收敛更快，且不容易陷入局部最优的情况。 此外，实验表明，学习率对神经网络的收敛也有比较大的影响。当学习率较大时，有可能收敛较快，当然也可能根本不收敛；当学习率较小时，容易进入局部最小值，且收敛较慢。还有就是网络深度不够或者神经元节点较少时，难以训练出较复杂的判决面(例如程序中想训练出一个圆形判决面，但是效果不好)。 如有精通矩阵的同学，可尝试推导一下损失函数对每个参数求偏导的矩阵求导写法。 接下来是“talk is cheap, show you the code”环节。 MATLAB代码及详细注释如下： 主函数：BP_xor_classifier.m clc;close all;clear; %%训练集 %% xor训练集 X = [0,0;0,1;1,0;1,1]'; %一列为一个样本 Y = [1,0;0,1;0,1;1,0]'; %标签一列为一个标签onehot编码后的 %%圆训练集（网络太浅，分类效果不好） % train_number = 400; % r_small = rand(1,train_number);%分布在小圆内的100个点 % r_big = r_small + rand(1,train_number); % angle_small = 2 * pi *rand(1,train_number); % angle_big = 2 * pi * rand(1,train_number); % small(1,:) = r_small .* cos(angle_small); % small(2,:) = r_small .* sin(angle_small); % small(3:4,:) = [ones(1,size(r_small,2));zeros(1,size(r_small,2))]; % big(1,:) = r_big .* cos(angle_big); % big(2,:) = r_big .* sin(angle_big); % big(3:4,:) = [zeros(1,size(r_small,2));ones(1,size(r_small,2))]; % joint_temp = [small,big]; % rerank_idx = randperm(size(joint_temp,2)); % joint = joint_temp(:,rerank_idx); % X = joint(1:2,:); % Y = joint(3:4,:); %%参考博客数据：测试用 % Y = [0.01;0.99]; % X = [0.05;0.1]; % theta_layer1 = [0.35,0.15,0.2;0.35,0.25,0.3]; % theta_layer2 = [0.6,0.4,0.45;0.6,0.5,0.55]; %%学习率和训练次数 learning_rate = 0.1; iter = 10000; %%搭建前向传播网络,初始化网络 layer1_unit = 2;%各层神经元的个数第一层为输入层，最后一层为输出层 layer2_unit = 2; layer3_unit = 2; %参数初始化 theta_layer1 = rand(layer2_unit,layer1_unit+1);%即M*N的随机矩阵，行数为下一层的神经元个数，一行为一组参数，即（b1,w1,w2） theta_layer2 = rand(layer3_unit,layer2_unit+1);%列数为前一层神经元个数+1（加上权重的参数1） for i = 1:iter for j = 1:size(X,2)%一列为一个样本，一个一个遍历该样本 %%搭建前馈式网络 %输入层 out_layer1 = [ones(1,size(X(:,j),2));X(:,j)];%往X的第一行插入一行1，即权重b的系数，用于后面的向量化相乘，这是第一层的输出（隐藏层的输入） %隐藏层 net_H = theta_layer1 * out_layer1; out_H = sigmoid(net_H); out_H_plus1 = [ones(1,size(out_H,2));out_H];%加一行,同上 %输出层 net_O = theta_layer2 * out_H_plus1; out_O = sigmoid(net_O); %%反向传播 %损失函数为MSE，即loss = (Y - out_O) .^ 2 / 2 %%首先求损失对theta_layer2求偏导（链式求导法则） %diff(loss/theta_layer2) = diff(loss/out_O) * diff(out_O/net_O) *diff(net_O/theta_layer2) %即损失对网络输出求导*网络输出对未激活的输出求导*未激活的输出对theta_layer2求导 %diff(loss/out_O) = -(Y - out_O ) %diff(out_O/net_O) = diff(sigmoid(net_O)) %diff(net_O/theta_layer2) = out_H_plus1 % theta_layer2_der = -(Y(:,j) - out_O) .* sigmoid_der(net_O) *out_H_plus1'; %上式中，diff(net_O/theta_layer2)（即out_H_plus1）中每一列的列向量， %分别是net_O对theta_layer2求的偏导结果 %这里做转置的原因是：前两项矩阵元素相乘的结果（2*1）， %每一项都要与out_H_plus1（3*1）中的元素相乘，得到一个2*3的矩阵 %%上式为矩阵写法，这里，可改写为以下程序更容易理解 % temp = -(Y(:,j) - out_O) .* sigmoid_der(net_O); % for m = 1:size(temp,1) % theta_layer2_der(m,:) = temp(m) * out_H_plus1; % end %%接下来对theta_layer1求偏导（同理，链式求导法则） %diff(loss/theta_layer1) = diff(loss/out_H) * diff(out_H/net_H) *diff(net_H/theta_layer1) %diff(loss/out_H) = diff(loss/out_O) * diff(out_O/net_O) *diff(net_O/out_H_plus1) %上式中，与diff(loss/theta_layer2)的区别只是最后一项是对out_H_plus1求的偏导 %diff(net_O/out_H_plus1) = theta_layer2'，矩阵乘法求导，Y = AX，则A需要转置 %上式中，diff(net_O/out_H_plus1)的结果包含偏置项的系数，在反向传播中不需要，将其舍去（结果的第一行） %diff(out_H/net_H) = diff(sigmoid(net_H)) %diff(net_H/theta_layer1) = out_layer1 % theta_layer1_der = theta_layer2(:,2:end)' * (-(Y(:,j) - out_O)) .*sigmoid_der(net_O) ... % .* sigmoid_der(net_H) *out_layer1'; %%另外的写法：记每一层的误差为Delta %diff(loss/net_O) = diff(loss/out_O) * diff(out_O/net_O) % Delta_3 = -(Y(:,j) - out_O) .* sigmoid_der(net_O); %损失函数为MSE时 Delta_3 = (out_O - Y(:,j))./((1 - out_O) .* out_O) .* sigmoid_der(net_O);%损失函数为交叉熵 %diff(loss/net_H) = diff(loss/out_O) * diff(out_O/net_O) *diff(net_O/out_H_plus1) * diff(out_H_plus1/net_H) %这里，需要将偏置那一项去掉，即diff(net_O/out_H_plus1) = theta_layer2(:,2:end)' Delta_2 = theta_layer2(:,2:end)' * Delta_3 .* sigmoid_der(net_H); theta_layer2_der = Delta_3 * out_H_plus1'; theta_layer1_der = Delta_2 * out_layer1'; %%更新参数 theta_layer2 = theta_layer2 - learning_rate * theta_layer2_der; theta_layer1 = theta_layer1 - learning_rate * theta_layer1_der; %注：这里两种方法计算所得的theta_layer1在小数点后第6位不相同，推测是由于软件精度所致 end %%计算损失并作图动态显示 % loss(i) = sum((Y(:,j) - out_O) .^ 2 / 2); loss(i) = sum(-Y(:,j) .* log(out_O) - (1 - Y(:,j)) .* log(1 - out_O)); % figure(1) % plot(i,loss,'<'); % drawnow; % hold on; end %%做预测 a = 0:0.01:1; b = 0:0.01:1; X = []; for i = 1:size(a,2) for j = 1:size(b,2) temp = [a(i);b(j)]; X = [X,temp]; end end % X = [0.7,10]'; for j = 1:size(X,2) out_layer1 =[ones(1,size(X(:,j),2));X(:,j)]; %隐藏层 net_H = theta_layer1 *out_layer1; out_H = sigmoid(net_H); out_H_plus1 =[ones(1,size(out_H,2));out_H];%加一行,同上 %输出层 net_O = theta_layer2 *out_H_plus1; predict(:,j) = sigmoid(net_O); end %%作图查看判决面（分类面） figure(1) plot(loss);%损失随迭代次数的变化 [~,c2]=max(predict); class_one_idx = find(c2 == 1); class_two_idx = find(c2 == 2); class_one = X(:,class_one_idx); class_two = X(:,class_two_idx); figure(2) plot(class_one(1,:),class_one(2,:),'ro'); hold on; plot(class_two(1,:),class_two(2,:),'b<'); sigmoid函数：sigmoid.m function a = sigmoid(X) for i = 1:size(X,1) for j = 1:size(X,2) a(i,j) = 1/(1+exp(-X(i,j)));%激活函数 end end end sigmoid的导数：sigmoid_der.m %%此函数为sigmoid的导数 function a = sigmoid_der(x) a = zeros(size(x)); for i = 1:size(x,1) a(i,:) = sigmoid(x(i,:)) *(eye(size(x,2)) - diag(sigmoid(x(i,:)))); end 参考： [1]https://mattmazur.com/2015/03/17/a-step-by-step-backpropagation-example/?spm=a2c4e.10696291.0.0.358f19a4xonKKs. [2]https://blog.csdn.net/zhaomengszu/article/details/77834845 [3]吴恩达机器学习视频
)
