《机器学习（周志华）》学习笔记（五）

Q：什么是人工神经网络？

人工神经网络是科学家模拟人类大脑的神经网络建立的数学模型。人工神经网络由一个个“人工神经元”组合而成。“人工神经元”也是一个数学模型，其本质是一个函数。所以人工神经网络的本质也是一个函数，而且是一个复杂的，包含很多变量和参数的函数。

Q：什么是人工神经元？

人工神经元是人工神经网络的基本单元，其本质也是一个函数。类似人类大脑的神经元，人工神经元模型也有树突、轴突、神经元中心等结构。最经典的神经元模型是“M-P神经元模型”。

一个神经元接收来自其他神经元传来的信号（变量），通过树突传输（参数）到神经元中心（参数），经过转换后（映射f），通过突触传递出去（函数值）。上图中红色的函数表达式就是一个人工神经元的数学表达。如果把上图的激活函数f设置成sigmoid函数，那么一个神经元本质上就是一个logistic 回归模型。

由上述神经元工作过程可知，一个神经元可以接受多个变量，输出一个结果。一个神经元的输出也可能是另一个神经元的输入。所以一个由多个神经元组成的神经网络，本质上就是一个复杂的、多层嵌套的复合函数。

Q:为什么神经网络要加入激活函数？

或问：神经网络引入非线性函数的用意是什么？

引入激活函数的根本目的是提高学习模型的学习（拟合数据）能力，其最初的动机是弥补线性模型学习能力不足的问题。

我们学习机器学习，都是从最简单的线性模型开始的（也有人从K临近算法开始）。我们知道线性模型虽然简单，但是功能却不差，起码预测一下波士顿的房价或者判断一下是否糖尿病是没问题的。但是线性模型说到底也只是一条直线或者一个直的平面，但是现实中的数据却不都是能用直线拟合或者能用平面分隔开的。比如依据培养时间预测细胞数目（函数图像为S型曲线，类似tanh）就不能用线性模型了，或者说线性模型的效果就很差了。

一句话，线性组合能力有限，不能表示复杂的规律；线性模型先天不足，不能学习复杂的情况。就像我们不能指望一个普通的6个月婴儿能掌握勾股定理，也不能指望一个普通的两岁小孩能理解量子力学。

所以我们在线性模型外面包上一层非线性函数，其实就是使其结构更加复杂，使得改造后的（广义）线性模型有能力表示复杂的映射，有能力学习到更复杂的情况，有能力拟合呈线性关系数据，也能拟合不呈线性关系的数据。

Q：神经网络有哪些分类？

按照神经元的层数来分，可以分为只有输入、输出层的单层网络（不计算输入层），代表是感知机；以及包含隐含层的多层网络。

按照网络结构来分，可以分为径向基函数网络（RBF）、竞争学习网络（ART）、自组织映射网络（SOM）、级联相关网络、递归神经网络（RNN）、Boltzman机等各种各样纷繁复杂的种类。

Q：感知机有什么作用？怎么训练一个感知机？

感知机是最简单的神经网络，只有输入层和输出层。从数学形式上看，感知机就是只有一个神经元的神经网络。一般来说，感知机用来进行二分类任务，或者实现逻辑“与”、“或”、“非”的操作。

感知机模型 $f(x)=sign(\vec{w}⋅\vec{x}+b)$ 的输出一般是二值的，即1或-1。

感知机的学习的学习规则很简单，也是类似于线性回归一样，使用梯度下降的思想，每读取一个样本，计算一次预测值，就调整一次参数（权值）。

$x_i$表示第i个输入，$w_i$表示第i个输入对应的参数（权值）。

Q：神经元如何组成神经网络？

通常来说，神经元会按层次结构组成一个网络，也就是一个多部图结构：每一层有任意数量个节点，但是同一层节点之间没有直接相连，只会和相邻层的结点直接相连，这也是经典的前馈神经网络的结构。

如此一来，本层的节点接收前一层节点的输出作为输入，经过计算后输出给后一层节点，每一层节点的输出结构都会是后一层节点的输入。

前馈神经网络只是最基本的网络架构，在此基础上有多种不同的变体：

竞争神经网络

竞争型网络的输出层之间会有直接联系，每个神经元都会和其他神经元竞争，胜利者才能顺利输出，落败者没有输出。典型的模型有自组织映射网络（SOM）。

级联相关网络

前馈网络的网络结构一般在训练开始前就设置好，训练过程中不会改变。级联相关网络则是可以在训练过程中动态改变网络结构：新增隐藏层和隐藏节点。

递归神经网络

前馈网络每一层的输入都是来自前一层的输出，而递归神经网络每一层的输入除了前一层的输出，还有本层在上一轮训练的输出，因此特别适合对序列数据（如文本数据）建模。Elman网络是常用的递归神经网络模型。

波兹曼机

波兹曼机只有两层，可见层和隐含层。标准的波兹曼机是一个全连接网络，每个神经元都和其他所有神经元相连，显然这种网络的复杂度就很高，难以训练。而受限波兹曼机则是同一层之间的神经元没有连接，类似一个二层前馈神经网络（不包含输入层）。

Q：神经网络训练过程是怎样的？

多层神经网络由于神经元层数比单层多，所以参数更多，也意味着学习能力更强，所以能够胜任图像处理、语音识别等任务。

最经典的多层网络训练算法是“误差逆传播算法（BP）”。其基本思想也是通过构建关于网络权重参数的损失函数来衡量神经网络的误差，然后通过梯度下降令损失函数取最小值，求得网络权重参数的最优值。

多层神经网络的参数也是按层划分的，所以参数的求解也要一层一层地进行。首先求出输出层的各个参数，求出倒数第二层的参数，再求解出倒数第三层的参数，以此类推，求解出所有层的参数。

由于参数是从输出层往后求解的，并且是通过损失函数计算出来的，因此成为误差反向传播算法。

下面是误差反向传播算法的数学推导

由于反向传播算法应用于复杂网络时会产生巨量的计算，为阐释清楚算法过程，这里把西瓜书上给出的示例网络进一步简化：整个网络只有输入层、一个隐含层，和输出层，每层只有一个神经元。训练样本也只有一个：$(x, y)$，其中x是一个实数（实际上应该是一维向量）。

在开始理解数学推导前，我们要先明确作者使用的数学符号所代表的含义：

• $x$ 训练集的样本数据。在这里只是一个实数。
• $y$ 训练集的标记数据，也就是实际结果。在这里只是一个实数（二分类任务的话就只有0或1）。
• $v$ 输入层和隐含层之间的连接权重（weights）
• $\gamma$  输入层和隐含层之间的连接偏置/阈值（bias)
• $\alpha = vx - \gamma$ 隐含层的线性组合计算值
• $\sigma$ Sigmoid函数：$\sigma(z) = \frac{1}{1 + e^{-z}}$
• $b$ 隐含层的输出结果，也是输出层的输入值
• $w$ 隐含层和输出层之间的连接权重（weights）
• $\theta$ 隐含层和输出层之间的连接偏置/阈值（bias)
• $\beta$ 输出层的线性组合计算值
• $\hat{y}$ 输出层的输出结果，也是神经网络的输出值，也就是预测结果 。

首先我们要拿到误差，然后才能反向传播

把训练样本$x$输入神经网络后，我们得到网络的预测值$\hat{y}$，然后我们可以用均方误差公式来计算本轮训练的误差：

$$E =  \frac{1}{2}  (y - \hat{y})^2$$

然后我们可以利用这个误差值来着手更新隐含层到输出层的连接权重。

更新隐含层到输出层的参数w和$\theta$

按照梯度下降算法，$w$和$\theta$的更新公式如下。其中$\eta$是学习率，一个我们自己定义的超参数。整个更新公式的含义就是让参数朝着学习误差E最小的方向（正是梯度的几何意义）移动一点，至于移动的步子多大，由学习率$\eta$决定。

$$w \leftarrow w + \eta\nabla w  = w + \eta \frac{\partial E}{\partial w}$$

$$\theta \leftarrow \theta + \eta\nabla \theta = \theta + \eta \frac{\partial E}{\partial \theta}$$

$\frac{\partial E}{\partial w}$和$\frac{\partial E}{\partial \theta}$这个又该怎么算？由于

$$E = \frac{1}{2} (y - \hat{y})^2 \\ \hat{y} = \sigma(\beta) \\ \beta = wb - \theta$$

通过导数的链式法则，有

$$\nabla w  = \frac{\partial E}{\partial w} = \frac{\partial E}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial \beta} \frac{\partial \beta}{\partial w}$$

$$\nabla \theta  = \frac{\partial E}{\partial \theta} = \frac{\partial E}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial \beta} \frac{\partial \beta}{\partial \theta}$$

其中

$$\frac{\partial E}{\partial \hat{y}} = \frac{\partial (\frac{1}{2} (y - \hat{y})^2 )}{\partial \hat{y}} =  y-\hat{y} \\ \frac{\partial \hat{y}}{\partial \beta}  = \frac{\partial (\frac{1}{1 + e^{-\beta}})}{\partial \beta} = \frac{1}{1 + e^{-\beta}}(1-\frac{1}{1 + e^{-\beta}}) = \hat{y}(1-\hat{y})\\ \frac{\partial \beta}{\partial w}  =  \frac{\partial ( bw - \theta)}{\partial w} =b$$

因此

$$\nabla w  = \frac{\partial E}{\partial w} = \frac{\partial E}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial \beta} \frac{\partial \beta}{\partial w} = (y-\hat{y}) · \hat{y}(1-\hat{y}) · b$$

$$\nabla \theta  = \frac{\partial E}{\partial \theta} = \frac{\partial E}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial \beta} \frac{\partial \beta}{\partial \theta} = (y-\hat{y}) · \hat{y}(1-\hat{y}) · 1$$

完成这一层后，我们接着更新后一层的参数，亦即更新输入层和隐含层之间的权重参数。

更新隐含层到输出层的参数v和$\gamma$

按照梯度下降算法，$v$和$\gamma$的更新公式如下。

$$v \leftarrow v + \eta\nabla v  = v + \eta \frac{\partial E}{\partial v}$$

$$\gamma \leftarrow \gamma + \eta\nabla \gamma = \gamma + \eta \frac{\partial E}{\partial \gamma}$$

$\frac{\partial E}{\partial v}$和$\frac{\partial E}{\partial \gamma}$又该怎么算？按照类似的流程，由于

$$E = (y - \hat{y})^2 \\ \hat{y} = \sigma(\beta) \\ \beta = wb - \theta \\ b = \sigma(\alpha) \\ \alpha = wx - \gamma$$

通过导数的链式法则，有

$$\nabla v  = \frac{\partial E}{\partial v} = \frac{\partial E}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial \beta} \frac{\partial \beta}{\partial b} \frac{\partial b}{\partial \alpha} \frac{\partial \alpha}{\partial v} = (y-\hat{y}) · \hat{y}(1-\hat{y}) · w · b(1-b) · x$$

$$\nabla \gamma  = \frac{\partial E}{\partial \gamma} = \frac{\partial E}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial \beta} \frac{\partial \beta}{\partial b} \frac{\partial b}{\partial \alpha} \frac{\partial \alpha}{\partial \gamma} = (y-\hat{y}) · \hat{y}(1-\hat{y}) · w · b(1-b) · 1$$

总结BP算法流程

对于上面我们给出的简化3层感知机，以及只有一个训练样本的情况，我们要循环往复地进行如下计算

$$\alpha = wx - \gamma \\ b = \sigma(\alpha) \\ \beta = wb - \theta \\ \hat{y} = \sigma(\beta) \\ w \leftarrow w + \eta(y-\hat{y})\hat{y}(1-\hat{y})b \\ \theta \leftarrow \theta + \eta(y-\hat{y})\hat{y}(1-\hat{y}) \\ v \leftarrow v + \eta(y-\hat{y})\hat{y}(1-\hat{y})wb(1-b)x  \\ \gamma \leftarrow \gamma + \eta(y-\hat{y})\hat{y}(1-\hat{y})wb(1-b)$$

一直计算到给定的循环次数，或者这一轮的误差与前一轮地误差相差小于一个阈值为止$E_t - E_{t-1} < \epsilon$。

以上便是误差反向传播算法在3层感知机，且每层只有一个神经元的情况下的数学推导。对于多层网络，每层多个神经元的情况，数学原理相同，具体可以参考吴恩达的机器学习课程或深度学习课程。

---

Talk is cheap, show me the code!

下面代码是按照西瓜书里的案例，实现的一个极其简陋的3层感知机分类器。这个模型没有momentum，没有penalty，没有adaptive，其他优化手段一切都没有，就只是一个简陋的网络结构。本来西瓜书里两层的激活函数都是sigmoid，后来测试的时候怎么测准确率都在0.56左右浮动，实在忍不了，把第一层的激活函数换了ReLU，这才好看一点，但仍然被其他模型吊打。可能是这份代码里面的bug词太多了吧？

"""
Naive 3 layer perceptron, uneffective, full of flaws!
 
:file: supervised.py
:author: Richy Zhu
:email: rickyzhu@foxmail.com
"""
import numpy as np
 
def binarize(y, threashold=0.5):
    '''
    Transform numeric data into binarize 0 and 1
 
    >>> y = np.array([0.27, 0.07, 0.56, 0.35, 0.32, 0.65])
    >>> binarize(y)
    array([0, 0, 1, 0, 0, 1])
    '''
    return np.where(np.array(y) < threashold, 0, 1)
 
def sigmoid(x): 
    '''
    Batch version Sigmoid function
 
    >>> x = np.array([-0.27, 0.07, 0.56, -0.35, 0.32, -0.65])
    >>> sigmoid(x)
    array([0.4329071 , 0.51749286, 0.63645254, 0.41338242, 0.57932425,
       0.34298954])
    '''
    return 1 / (1 + np.exp(-x))
 
def relu(x):
    '''
    Batch version ReLU function
 
    >>> x = np.array([-0.27, 0.07, 0.56, -0.35, 0.32, -0.65])
    >>> relu(x)
    array([0.  , 0.07, 0.56, 0.  , 0.32, 0.  ])
    '''
    return np.where(x<0, 0, x)
 
 
class My3LPClassifier:
    '''An naive 3 layer perceptron classifier'''
    
    def __init__(self, layer_dims):
 
        assert len(layer_dims)==3, 'layer_dims should be 3 elements'
        self.input_size = layer_dims[0]
        self.hidden_size = layer_dims[1]
        self.output_size = layer_dims[-1]
        
        self.V = np.random.randn(self.input_size, self.hidden_size)
        self.gamma = np.zeros([1, self.hidden_size])
        self.W = np.random.randn(self.hidden_size, self.output_size)
        self.theta = np.zeros([1, self.output_size])
        
    def _forward(self, x):
        '''forward computing'''
        x = x.reshape([1,-1])
        alpha = np.dot(x, self.V) - self.gamma
        b =     relu(alpha)
        # b =     sigmoid(alpha)
        beta =  np.dot(b, self.W)
        y_hat = sigmoid(beta) - self.theta
        return alpha, b, beta, y_hat
    
    def _loss(self, y, y_hat):
        '''mean square error'''
        return np.sum(np.square(y-y_hat)) / len(y)
 
    def _back(self, x, y, alpha, b, beta, y_hat, eta):
        '''back propagation'''
        g = y_hat * (1-y_hat) * (y - y_hat)
        self.W = self.W + eta * (g * b.reshape([-1,1]))
        self.theta = self.theta - eta * g
        e = np.where(alpha<0, 0, alpha) * np.sum(self.W * g, axis=1)
        # e = b * (1-b) * np.sum(self.W * g, axis=1)
        self.V = self.V + eta * (e * x.reshape([-1,1]))
        self.gamma = self.gamma - eta * e
 
    def fit(self, X, y, eta=1e-5, eps=1e-6, max_iter=10000):
        '''
        Build an random forest classifier
 
        Parameters
        ----------
        X: ndarray of shape (m, n)
            sample data where row represent sample and column represent feature
        y: ndarray of shape (m, k)
            labels of sample data
 
        Returns
        -------
        self
            trained model
        '''
        self.X = X
        self.y = y
        
        i = 0
        prev_loss = np.inf
        while i < max_iter:
            predicted = []
            for x, y in zip(self.X, self.y):
                alpha, b, beta, y_hat = self._forward(x)
                self._back(x, y, alpha, b, beta, y_hat, eta)
                predicted.append(y_hat)
            this_loss = self._loss(self.y, np.array(predicted))
            if prev_loss - this_loss < eps:
                break
            prev_loss = this_loss
            i += 1
        return self
    
    def predict(self, X):
        '''
        Make prediction by the trained model.
 
        Parameters
        ----------
        X: ndarray of shape (m, n)
            data to be predicted, the same shape as trainning data
 
        Returns
        -------
        C: ndarray of shape (m, k)
            Predicted class label(s) per sample.
        '''
        y = []
        for x in X:
            y.append(self._forward(x)[-1])
        return binarize(np.array(y))
    

测试一下

import numpy as np
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.metrics import  accuracy_score
 
print('\nMultilayer Perceptron')
print('---------------------------------------------------------------------')
X, y = make_classification(n_samples=1000, n_features=4,
                           n_informative=2, n_redundant=0,
                           random_state=0, shuffle=False)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
mymlp = My3LPClassifier(layer_dims=[X_train.shape[1], 3, 1])
mymlp.fit(X_train, y_train)
print('My 3LP:',  accuracy_score(mymlp.predict(X_test).flatten(), y_test))
 
from sklearn.neural_network import MLPClassifier
skmlp = MLPClassifier()
skmlp.fit(X_train, y_train)
print('Sk 3LP:',  accuracy_score(skmlp.predict(X_test), y_test))

测试结果如下，可谓惨不忍睹，scikit-learn里面任何一个分类器拿出来都能吊打，侧面反映出网络的优化是多么重要。

Multilayer Perceptron
---------------------------------------------------------------------
My 3LP: 0.76
Sk 3LP: 0.96

更多代码请参考https://github.com/qige96/programming-practice/tree/master/machine-learning

---

本作品首发于简书 和 博客园平台，采用知识共享署名 4.0 国际许可协议进行许可。