多层感知机MLP多分类python实践

2019-04-12

之前不太明白神经网络什么层和反向传播什么的，这两天试着写了个纯python+numpy版本的简单神经网络(也就是多层感知机)，感觉对原理清楚了不少。
首先生成一些不能线性分割的数据，分三类以上，使用softmax作为输出函数来进行分类：

import numpy as np
import matplotlib.pyplot as plt
N = 100 # 每一个类别的生成的点的数量
D = 2 # 每个点的维度，这里使用平面，所以是2维数据
K = 3 # 类别数量，我们一共生成3个类别的点

# 所有的样本数据，一共300个点，每个点用2个维度表示
# 所有训练数据就是一个300*2的二维矩阵
X = np.zeros((N*K, D))
# 标签数据，一共是300个点，每个点对应一个类别，
# 所以标签是一个300*1的矩阵
y = np.zeros(N*K, dtype='uint8')
np.random.seed(0)
# 生成训练数据
for j in range(K):
    ix = range(N*j, N*(j+1))
    r = np.linspace(0.0, 1, N)
    t = np.linspace(j*4, (j+1)*4, N) + np.random.randn(N)*0.2
    X[ix] = np.c_[r*np.sin(t), r*np.cos(t)]
    y[ix] = j
    
plt.scatter(X[:, 0], X[:, 1], c=y, s=40, cmap=plt.cm.Spectral)
plt.show()

打印出来是这样的

接着定义一些会用到的函数，包括激活函数和它的导数，以及打印模型的函数等等：

def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))

def d_sigmoid(y):
    return y * (1.0 - y)

def tanh(x):
    return np.tanh(x)

def d_tanh(y):
    return 1. - y**2


def softmax(Z):
        return np.exp(Z) / np.sum(np.exp(Z), axis=1, keepdims=True)
    
def cross_entropy(yhat,y):
    return -y*np.log(yhat)-(1-y)*np.log(1-yhat)

# 一个函数来画决策边界
def plot_decision_boundary(data,label,pred_func):
 
    # 设定最大最小值，附加一点点边缘填充
    x_min, x_max = data[:, 0].min() - .5, data[:, 0].max() + .5
    y_min, y_max = data[:, 1].min() - .5, data[:, 1].max() + .5
    h = 0.01
 
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))
 
    # 用预测函数预测一下
    Z = pred_func(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)
 
    # 然后画出图
    plt.contourf(xx, yy, Z, cmap=plt.cm.Spectral)
    plt.scatter(data[:, 0], data[:, 1], c=label, cmap=plt.cm.Spectral)

先看下只用softmax函数(相当于没有隐藏层的神经网络)做分类的效果：

看到只有直线作为决策边界，正确率也只有51%
接下来就是使用神经网络的训练过程了，先来一些原理。简单来说就是每一轮都先计算出一个分类结果(前向传播)，然后与标准答案做对比计算出损失，再由损失反过来修改参数的权重(反向传播)。

首先看前向传播的过程，这个很简单，跟逻辑回归类似。一共有三层，输入层，隐藏层和输出层，输入层传进去的数据X和权重W做矩阵乘法得到Z，之后Z通过激活函数f(sigmoid/tanh/relu)得到A，这样就有了隐藏层的输出，再把隐藏层的输出作为输出层的输入，同样的过程得到Z后通过softmax函数(不是激活函数了)得到预测的结果yhat。代码也很简单，这里考虑了不同的激活函数

#forward
    Z1 = np.dot(X,self.W1) + self.b1
    if self.activation == 'sigmoid':
        A1 = sigmoid(Z1)
    elif self.activation == 'tanh':
        A1 = tanh(Z1)
    else:
        #relu
        A1 = np.maximum(0,Z1)            
    Z2 = np.dot(A1,self.W2) + self.b2
    yhat = softmax(Z2)

接着是反向传播的过程，首先要计算出损失，这样才可以对损失求梯度来减小损失。softmax都使用交叉熵作为损失函数，因为这样求导后很简单。具体过程可以看这篇文章https://zhuanlan.zhihu.com/p/37740860
得到了输出层的导数是yhat-y，之后由链式法则(这写的参数和上边那个图不一样，A2就是yhat)

1	dJ/dW2 = dJ/dA2 * dA2/dZ2 * dZ2/dW2

前两项就是输出层的导数delta3，第三项就是A1
继续展开

1
2
3

dJ/dW1 = dJ/dA1 * dA1/dZ1 * dZ1/dW1
	   = dJ/dA2 * dA2/dZ2 * dZ2/dA1 * dA1/dZ1 * dZ1/dW1
	   = delta3 * W2 * d_active(Z1) * X

dZ2/dA1就是W2，dA1/dZ1就是激活函数的导数，dZ1/dW1就是X
和这张图对照着看，Z和A的下标需要减1

这样代码也很好写了

#loss
    loss = np.mean(cross_entropy(yhat,y))
#   reg_loss = 0.5*reg*np.sum(self.W1*self.W1) + 0.5*reg*np.sum(self.W2*self.W2)
#   loss += reg_loss
#backward
	delta3 = yhat - y
#	y不进行onehot化的话按下面这两行这么写，等价于yhat-y
#   delta3 = yhat
#   delta3[range(X.shape[0]),y] -= 1
    delta3 /= X.shape[0] #学习率较大的时候必须加，否则不收敛

    dJdW2 = np.dot(A1.T,delta3)
    dJdb2 = np.sum(delta3, axis=0, keepdims=True)
    if self.activation == 'sigmoid':
        delta2 = np.dot(delta3,self.W2.T)*d_sigmoid(A1)
    elif self.activation == 'tanh':
        delta2 = np.dot(delta3,self.W2.T)*d_tanh(A1)
    else:
        #relu 导数
        delta2 = np.dot(delta3,self.W2.T)
        delta2[A1 <= 0] = 0

    dJdW1 = np.dot(X.T,delta2)
    dJdb1 = np.sum(delta2, axis=0)
#   dJdW2 += reg * dJdW2
#   dJdW1 += reg * dJdW1
    self.W1 -= self.learning_rate * dJdW1
    self.b1 -= self.learning_rate * dJdb1
    self.W2 -= self.learning_rate * dJdW2
    self.b2 -= self.learning_rate * dJdb2

这样训练函数就写好了，最后预测对softmax的结果取argmax就可以了。完整代码：

def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))

def d_sigmoid(y):
    return y * (1.0 - y)

def tanh(x):
    return np.tanh(x)

def d_tanh(y):
    return 1. - y**2


def softmax(Z):
        return np.exp(Z) / np.sum(np.exp(Z), axis=1, keepdims=True)
    
def cross_entropy(yhat,y):
    return -y*np.log(yhat)-(1-y)*np.log(1-yhat)

class MLPClassifier():
    """多层感知机，BP 算法训练"""

    def __init__(self,hid_neurals=10,epochs=1000, batch_size=1, learning_rate=0.01, activation='relu'):
        """
        :param layers: 网络层结构
        :param activation: 激活函数
        :param epochs: 迭代轮次
        :param learning_rate: 学习率 
        """
        self.epochs = epochs
        self.learning_rate = learning_rate
        self.hid_neurals = hid_neurals
        self.batch_size = batch_size
        self.activation = activation
#         self.reg = 1e-3 # regularization strength
        
    def fit(self,X,y):
        self.W1 = np.random.randn(X.shape[1], self.hid_neurals)
        self.b1 = np.zeros((1, self.hid_neurals))
        self.W2 =  np.random.randn(self.hid_neurals,y.shape[1])#分类个数
        self.b2 = np.zeros((1, y.shape[1]))
        for i in range(self.epochs):
            #forward
            Z1 = np.dot(X,self.W1) + self.b1
            if self.activation == 'sigmoid':
                A1 = sigmoid(Z1)
            elif self.activation == 'tanh':
                A1 = tanh(Z1)
            else:
                #relu
                A1 = np.maximum(0,Z1)            
            Z2 = np.dot(A1,self.W2) + self.b2
            yhat = softmax(Z2)
 
            #loss
            loss = np.mean(cross_entropy(yhat,y))
#             reg_loss = 0.5*reg*np.sum(self.W1*self.W1) + 0.5*reg*np.sum(self.W2*self.W2)
#             loss += reg_loss
            
            if i % 1000 == 0:
                print ("iteration %d loss %f" % (i, loss))
            
            #backward
            delta3 = yhat - y
#             print(yhat[0])
            
            delta3 /= X.shape[0] #学习率较大的时候必须加，否则不收敛
#             print(delta3[0])
            
            dJdW2 = np.dot(A1.T,delta3)
            dJdb2 = np.sum(delta3, axis=0, keepdims=True)
            if self.activation == 'sigmoid':
                delta2 = np.dot(delta3,self.W2.T)*d_sigmoid(A1)
            elif self.activation == 'tanh':
                delta2 = np.dot(delta3,self.W2.T)*d_tanh(A1)
            else:
                #relu 导数
                delta2 = np.dot(delta3,self.W2.T)
                delta2[A1 <= 0] = 0
            
            dJdW1 = np.dot(X.T,delta2)
            dJdb1 = np.sum(delta2, axis=0)
#             dJdW2 += reg * dJdW2
#             dJdW1 += reg * dJdW1
            self.W1 -= self.learning_rate * dJdW1
            self.b1 -= self.learning_rate * dJdb1
            self.W2 -= self.learning_rate * dJdW2
            self.b2 -= self.learning_rate * dJdb2
    
    def predict(self,X):
        Z1 = np.dot(X,self.W1) + self.b1
        if self.activation == 'sigmoid':
            A1 = sigmoid(Z1)
        elif self.activation == 'tanh':
            A1 = tanh(Z1)
        else:
            A1 = np.maximum(0,Z1)
        Z2 = np.dot(A1,self.W2) + self.b2
        yhat = softmax(Z2)
        return np.argmax(yhat, axis=1)

from sklearn import preprocessing
mlp =  MLPClassifier(epochs=10000,learning_rate=0.1,activation='tanh')
enc = preprocessing.OneHotEncoder()
target_train = enc.fit_transform(y.reshape(-1, 1)).toarray()  # 对目标集独热编码
mlp.fit(X,target_train)
plot_decision_boundary(X,y,lambda x:mlp.predict(x))
np.mean(mlp.predict(X)==y)

正确率97.7%，还是很有效果的