神经网络笔记(Neural Network)

给定训练样本集(x(i),y(i))，神经网络定义了一种参数为W,b的复杂非线性的假设hW,b(x)，来拟合数据。

最简单的神经网络是神经元，如下图：

这个神经元输入x1,x2,x3和一个+1截距(b表示)，输出为hW,b(x)=f(WTx)=f(∑3i=1Wixi+b)。函数f:R↦R称为激活函数。常用的激活函数为sigmoid函数：f(z)=11+e−z，f′(z)=f(z)(1−f(z)，而如此，单个神经元就直接是一个logistic regression。
其他常用的激活函数还有tanh函数：f(z)=tanh(z)=ez−e−zez+e−z，f′(z)=1−(f(z))2，它将sigmoid函数的值域扩展到了[−1,1]；以及最近研究发现的rectified linear函数：f(z)=max(0,z)，其导数当z≤0时为0，其余为1。
以下是三个函数的图像：

Neural Network model

将多个神经元连接起来便成了神经网络，如下图的例子：

+1称为偏置节点，截距项。最左边一层为输入层，最右层为输出层。中间层为隐藏层，因为不能在训练样本集中观测到他们的值。不包括偏置结点，以上神经网络拥有3个输入单元，3个隐藏单元和1个输出单元。

用nl来表示神经网络的层数，比如以上例子nl=3；用Ll来表示第l层，比如以上例子输入层为L1，Lnl为输出层。
以上例子中，神经网络的参数为(W,b)=(W(1),b(1),W(2),b(2))，W(l)ij表示第l层第j单元与第l+1层第i单元的偏置项。在以上神经网络中W(1)∈R3∗3,W(2)∈R1∗3。用sl表示第l层的节点数，不计偏置单元。

用a(l)i表示第l层第i单元的激活值。当l=1时，a(1)i=xi，为第i个输入值。对于给定参数集合W,b，神经网络就可以按照函数hW,b(x)来计算输出结果。

Forward Propagation

以上神经网络计算步骤如下：
a(2)1=f(W(1)11x1+W(1)12x2+W(1)13x3+b(1)1)
a(2)2=f(W(1)21x1+W(1)22x2+W(1)23x3+b(1)2)
a(2)3=f(W(1)31x1+W(1)32x2+W(1)33x3+b(1)3)
hW,b(x)=a(3)1=f(W(2)11a(2)1+W(2)12a(2)2+W(2)13a(2)3+b(2)1)

在上面等式中，通常用z(l)i表示l层第i单元总的带权输出。比如z(2)i=∑nj=1W(1)ijxj+b(1)i，因此a(l)i=f(z(l)i)。

如果将函数f()扩展到向量，比如f([z1,z2,z3])=[f(z1),f(z2),f(z3)]，那么，以上等式可以表达得更简洁：
z(2)=W(1)x+b(1)
a(2)=f(z(2))
z(3)=W(2)a(2)+b(2)
hW,b(x)=a(3)=f(z(3))

这个步骤叫做前向传播(Forward propagation)。
更一般地，定义a(1)=x表示输入层的值，则从l层向l+1的前向传播过程：
z(l+1)=W(l)a(l)+b(l)
a(l+1)=f(z(l+1))

Backpropagation Alrithm

假设我们有训练集{(x(1),y(1)),...,(x(m),y(m))}，我们可以使用batch gradient descent来训练神经网络。
具体来说，对于单个训练样例(x,y)，cost function定义为J(W,b;x,y)=12||hW,b(x)−y||2。
对于m个样例的数据集，cost function定义为：
J(W,b)=[1m∑mi=1J(W,b;x(i),y(i))]+λ2∑nl−1l=1∑sli=1∑sl+1j=1(W(l)ji)2=[1m∑mi=1(12||hW,b(x(i))−y(i)||2)]+λ2∑nl−1l=1∑sli=1∑sl+1j=1(W(l)ji)2

第一项为均方差项，第二项为规则化项，减小权重幅度，防止过拟合。

对于参数W和b，应该用一个很小的，接近零的随机值，比如正态分布Normal(0,ϵ2)生成随机值，其中ϵ设置为0.01。

因为J(W,b)是一个非凸函数，梯度下降法很可能会收敛到局部最优解；但是在实际应用中，梯度下降法通常能得到令人满意的结果。

最后，需要再次强调的是，要将参数进行随机初始化，而不是全部置为0。如果所有参数都用相同的值作为初始值，那么所有隐藏层单元最终会得到与输入值有关的、相同的函数（也就是说，对于所有i，W(1)ij都会取相同的值，那么对于任何输入x都会有：a(2)1=a(2)2=a(2)3=…）。随机初始化的目的是使对称失效。

梯度下降法中每一次迭代都按照如下公式对参数 \textstyle W 和\textstyle b 进行更新：
W(l)ij=W(l)ij−α∂∂W(l)ijJ(W,b)
b(l)i=b(l)i−α∂∂b(l)iJ(W,b)

其中α是学习速率。其中关键步骤是计算偏导数。

接着介绍反向传播算法，它是一种计算偏导的有效方法。

用反向传播算法来计算 ∂∂W(l)ijJ(W,b;x,y)和∂∂b(l)iJ(W,b;x,y)，这两项是单个样例(x,y)的代价函数J(W,b;x,y)的偏导数。求出偏导数，就可以推导出整体代价函数J(W,b)的偏导数：
∂∂W(l)ijJ(W,b)=[1m∑mi=1∂∂W(l)ijJ(W,b;x(i),y(i))]+λW(l)ij
∂∂b(l)iJ(W,b)=1m∑mi=1∂∂b(l)iJ(W,b;x(i),y(i))

以上两行公式稍有不同，第一行比第二行多出一项，是因为权重衰减是作用于W而不是b。

反向传播算法主要思想是先进行向前传播，计算网络中所有的值。然后针对第l层的每一个节点i，计算出其“残差”δ(l)i，该残差表明了该节点对最终输出值的残差产生了多少影响。对于最终的输出节点，可以直接算出网络产生的激活值与实际值之间的差距，将这个差距定义为δ(nl)i（第nl 层表示输出层）。

具体实现：

1. 进行前馈传导计算，利用前向传导公式，得到L2,L3,…直到输出层Lnl的激活值。
2. 对于第nl层（输出层）的每个输出单元i，我们根据以下公式计算残差：
   δ(nl)i=∂∂znliJ(W,b;x,y)=∂∂znli12||y−hW,b(x)||2=∂∂znli12∑Snlj=1(yj−a(nl)j)2=∂∂znli12∑Snlj=1(yj−f(z(nl)j))2=−(yi−f(z(nl)i))⋅f′(z(nl)i)=−(yi−a(nl)i)⋅f′(z(nl)i)
3. 对l=nl−1,nl−2,nl−3,…,2的各个层，第l层的第i个节点的残差计算方法如下：δ(nl−1)i=∂∂znl−1iJ(W,b;x,y)=∂∂znl−1i12∥∥y−hW,b(x)∥∥2=∂∂znl−1i12∑Snlj=1(yj−a(nl)j)2=12∑Snlj=1∂∂znl−1i(yj−a(nl)j)2=12∑Snlj=1∂∂znl−1i(yj−f(z(nl)j))2=∑Snlj=1−(yj−f(z(nl)j))⋅∂∂z(nl−1)if(z(nl)j)=∑Snlj=1−(yj−f(z(nl)j))⋅f′(z(nl)j)⋅∂z(nl)j∂z(nl−1)i=∑Snlj=1δ(nl)j⋅∂z(nl)j∂znl−1i=∑Snlj=1(δ(nl)j⋅∂∂znl−1i∑Snl−1k=1f(znl−1k)⋅Wnl−1jk)=∑Snlj=1δ(nl)j⋅Wnl−1ji⋅f′(znl−1i)=(∑Snlj=1Wnl−1jiδ(nl)j)f′(znl−1i)
   将上式中的nl−1与nl的关系替换为l与l+1的关系，就可以得到：δ(l)i=(∑sl+1j=1W(l)jiδ(l+1)j)f′(z(l)i)
   4.计算偏导数：∂∂W(l)ijJ(W,b;x,y)=a(l)jδ(l+1)i∂∂b(l)iJ(W,b;x,y)=δ(l+1)i

扩展函数f(x)与f′(x)到向量，即f′([z1,z2,z3])=[f′(z1),f′(z2),f′(z3)]，可以将上面的步骤简化表示如下：

1. 进行前馈传导计算，利用前向传导公式，得到L2,L3,…直到输出层Lnl的激活值。
2. 对输出层（第nl层），计算：δ(nl)=−(y−a(nl)).∗f′(z(nl))。
3. 对于l=nl−1,nl−2,nl−3,…,2的各层，计算：δ(l)=((W(l))Tδ(l+1)).∗f′(z(l))。
4. 计算最终需要的偏导数值：
   ∇W(l)J(W,b;x,y)=δ(l+1)(a(l))T
   ∇b(l)J(W,b;x,y)=δ(l+1)。

Vectorization

使用学习算法时，一段更快的代码通常意味着项目进展更快。
总而言之，矢量化编程是提高算法速度的一种有效方法。思想是尽量使用被高度优化的数值运算操作来实现学习算法。
在matlab中，矢量化的诀窍在于：代码中尽量避免显式的for循环。

调试方法：刚开始编写程序的时候，你可能会选择不使用太多矢量化技巧来实现你的算法，并验证它是否正确（可能只在一个小问题上验证）。在确定它正确后，你可以每次只矢量化一小段代码，并在这段代码之后暂停，以验证矢量化后的代码计算结果和之前是否相同。

Logistic Regression

用批量梯度上升法对logistic回归分析模型进行训练，其模型如下：
hθ(x)=11+exp(−θTx)
符号定义：设x0=1，于是x∈Rn+1，θ∈Rn+1，θ0为截距。假设我们有m个训练样本{(x(1),y(1)),...,(x(m),y(m))}，而批量梯度上升法的更新法则是：θ:=θ+α∇θl(θ),这里的l(θ)是对数似然函数，∇θl(θ)是其导函数。而∇θℓ(θ)=∑mi=1(y(i)−hθ(x(i)))x(i)j。

x(:,i)代表第i个训练样本x(i)，x(j,i)就代表x(i)j，表示第i个训练样本向量的第j个元素。
y表示由训练样本集合的全体类别标号所构成的行向量，则该向量的第i个元素y(i)就代表上式中的y(i)∈{0,1}。
初实现，循环嵌套，速度非常慢：

grad = zeros(n+1,1);for i = 1 : m,  h = sigmoid(theta'*x(:,i));  temp = y(i) - h;   for j = 1 : n+1,    grad(j) = grad(j) + temp * x(j,i);   end;end;

优化内层循环：

grad = zeros(n+1,1);for i=1:m,  grad = grad + (y(i) - sigmoid(theta'*x(:,i)))* x(:,i);end;

彻底优化：
特别的，假定b是一个列向量，A是一个矩阵，我们用以下两种方式来计算A*b：

% 矩阵-向量乘法运算的低效代码grad = zeros(n+1,1);for i=1:m,  grad = grad + b(i) * A(:,i);  % 通常写法为A(:,i)*b(i)end;% 矩阵-向量乘法运算的高效代码grad = A*b;

将b(i)看成(y(i) - sigmoid(theta’*x(:,i)))，A看成x，我们就可以使用以下高效率的代码：

grad = x * (y- sigmoid(theta'*x));

这里假定sigmoid(z)函数接受一个向量形式的输入z，依次对输入向量的每个元素施行sigmoid函数，最后返回运算结果，因此sigmoid(z)的输出结果是一个与z有相同维度的向量。

Neural Network

Forward Propagation

考虑一个三层网络(一个输入层、一个隐含层、以及一个输出层)，并且假定x是包含一个单一训练样本x^{(i)} \in \Re^{n} 的列向量。则向量化的正向传播步骤如下：
z(2)=W(1)x+b(1)a(2)=f(z(2))z(3)=W(2)a(2)+b(2)hW,b(x)=a(3)=f(z(3))
这对于单一训练样本而言是非常有效的一种实现，但是当我们需要处理m个训练样本时，则需要把如上步骤放入一个for循环中。
x表示包含输入训练样本的矩阵，x(:,i)代表第i个训练样本。则x正向传播步骤可如下实现：

% 非向量化实现for i=1:m,   z2 = W1 * x(:,i) + b1;  a2 = f(z2);  z3 = W2 * a2 + b2;  h(:,i) = f(z3);end;

向量化：

% 正向传播的向量化实现z2 = W1 * x + repmat(b1,1,m);a2 = f(z2);z3 = W2 * a2 + repmat(b2,1,m);h = f(z3)

其中repmat(A,m,n)函数能够把矩阵A，扩展成m行n列个A。比如：
B=repmat( [1 2;3 4],2,3)得到：
B=⎛⎝⎜⎜⎜131324241313242413132424⎞⎠⎟⎟⎟

激活函数实现：

% 低效的、非向量化的激活函数实现function output = unvectorized_f(z)output = zeros(size(z))for i=1:size(z,1),   for j=1:size(z,2),    output(i,j) = 1/(1+exp(-z(i,j)));  end; end;end% 高效的、向量化激活函数实现function output = vectorized_f(z)output = 1./(1+exp(-z));     % "./" 在Matlab或Octave中表示对矩阵的每个元素分别进行除法操作end

Back Propagation

假定网络的输出有s3维，因而每个样本的类别标号向量就记为y(i)∈Rs3，其中第i列y(:,i)就是y(i)。
还是以三层为例子，优化前：

% 'for' version   for i = 1 : m       y_vec = zeros(1, num_labels);       y_vec(y(i)) = 1;       delta3 = a3(i, :) - y_vec; %1 10       delta2 = delta3 * Theta2 .* [0 sigmoidGradient(z2(i, :))]; %1 26       delta2 = delta2(2 : end); %1 25       Theta2_grad = Theta2_grad + delta3' * a2(i, :);       Theta1_grad = Theta1_grad + delta2' * a1(i, :);   end

向量化：

% vectorize    y_vec = zeros(m, num_labels);    for i = 1 : m        y_vec(i, y(i)) = 1;    end    delta3 = a3 - y_vec;    delta2 = delta3 * Theta2 .* [zeros(m, 1), sigmoidGradient(z2)];    delta2 = delta2(:, 2:end);    Theta2_grad = delta3' * a2;    Theta1_grad = delta2' * a1;

summary

拿之前那个手写识别来试了下，这个优化是相当可观的。

向量化之前：

向量化之后：

优化之后的costFunction：

function [J grad] = CostFunction(nn_params, ...                                   input_layer_size, ...                                   hidden_layer_size, ...                                   num_labels, ...                                   X, y, lambda)    m = size(X, 1);    Theta1 = reshape(nn_params(1:hidden_layer_size * (input_layer_size + 1)), ...                     hidden_layer_size, (input_layer_size + 1));    Theta2 = reshape(nn_params((1 + (hidden_layer_size * (input_layer_size + 1))):end), ...                     num_labels, (hidden_layer_size + 1));    J = 0;    Theta1_grad = zeros(size(Theta1));    Theta2_grad = zeros(size(Theta2));% =========== Feedforward ========== %    a1 = [ones(m, 1) X]; %5000 401    z2 = a1 * Theta1';   %5000 25    a2 = [ones(m, 1) sigmoid( z2 )]; %5000 26    z3 = a2 * Theta2';               %5000 10    a3 = sigmoid( z3 );              %5000 10    hx = a3;                         %5000 10% ====================================== %% ===========  cost =================== %% 'for' version%    for i = 1 : m%        y_vec = zeros(1, num_labels);%        y_vec(y(i)) = 1;%        J = J + sum( -y_vec .* log(hx(i, :)) - (1 - y_vec) .* log(1 - hx(i, :)) );%    end% vectorize   y_vec = zeros(m, num_labels);   for i = 1 : m       y_vec(i, y(i)) = 1;   end   tmp = -y_vec .* log(hx) - (1 - y_vec) .* log(1 - hx);   J = J + sum( tmp(:) );    J = J / m;    J = J + lambda/(2*m) * (sum(sum(Theta1(:,2:end).^2))+sum(sum(Theta2(:,2:end).^2)));% ======================================= % % =========== backpropagation ========== %% 'for' version%    for i = 1 : m%        y_vec = zeros(1, num_labels);%        y_vec(y(i)) = 1;%        delta3 = a3(i, :) - y_vec; %1 10%        delta2 = delta3 * Theta2 .* [0 sigmoidGradient(z2(i, :))]; %1 26%        delta2 = delta2(2 : end); %1 25%        Theta2_grad = Theta2_grad + delta3' * a2(i, :);%        Theta1_grad = Theta1_grad + delta2' * a1(i, :);%    end% vectorize    y_vec = zeros(m, num_labels);    for i = 1 : m        y_vec(i, y(i)) = 1;    end    delta3 = a3 - y_vec;    delta2 = delta3 * Theta2 .* [zeros(m, 1), sigmoidGradient(z2)];    delta2 = delta2(:, 2:end);    Theta2_grad = delta3' * a2;    Theta1_grad = delta2' * a1;% ====================================== %    Theta2_grad = Theta2_grad / m;    Theta1_grad = Theta1_grad / m;% regularization    Theta1(:, 1) = 0;    Theta1_grad = Theta1_grad + lambda / m * Theta1;    Theta2(:, 1) = 0;    Theta2_grad = Theta2_grad + lambda / m * Theta2;    grad = [Theta1_grad(:) ; Theta2_grad(:)];end