《统计学习方法》之感知机模型

最近在重新看李航的统计学习方法，总结下每章的内容，并使用python复现。

基本概念

感知机定义 ：

输入空间 $\chi \subseteq R^n$ ，输出空间是 $Y=\{+1,-1\}$ ，输入 $x \in \chi$ ，表示实例的特征向量，对应于输入空间的点；输出 $y \in Y$ ，是该实例的类别。有输入空间到输出空间的映射函数决定：
$f(x)=sign(w \cdot x+b)$
其中 $sign(x)$ 表示符号函数。 $x$ 非负时为+1，否则是-1。

注：维基百科上 $sign(x)$ 函数在 $x=0$ 处取值定义为0，该函数简写为 $sgn(x)$ ，但《统计学习方法》中符号函数在0处取值定义为1，《机器学习》（周志华）中p98定义 $sgn(x)$ 为阶跃函数，在所有自变量非正时都取0。

感知机是线性分类模型，判别模型。可以用超平面 $wx+b=0$ 解释，该平面将空间分为两部分，即正负两类别。

数据集的线性可分性：

是指存在某个超平面 $wx+b=0$ 可以完整地将数据集的正负实例点划分到两边。对 $y_i=+1$ 的样本，都有 $w_ix+b>0$ ，而所有 $y_i=-1$ 的样本，都有 $w_ix+b<0$ 。这种数据集称为线性可分数据集。

感知机的学习策略：

输入空间中的任意点 $x_0$ 到超平面 $S$ 的距离定义是 $\frac 1 {||w||} |wx_0+b|$ ，其中 $||w||$ 表示 $w$ 的 $L_2$ 范数。

所有误分类的数据都满足不等式 $-y_i(wx_i+b)>0$ ，这是显然的，因为实例类别和计算的值异号。而 $y_i$ 的取值是 $\pm1$ ，因此误分类点到超平面的距离可以写成 $-\frac 1 {||w||} y_i(wx_i+b)$ ，距离公式里绝对值是正时， $-y_i$ 是正1，反之是负1，从而把绝对值符号去掉。

数据集中所有的误分类点集合是 $M$ ，则全部误分类点到超平面的总距离是
$-\frac 1 {||w||} \sum \limits_{x_i \in M}y_i(wx_i+b)$
不考虑前面的 $||w||$ ，定义感知机的损失函数为 $L(w,b)= -\sum \limits_{x_i \in M}y_i(wx_i+b)$

损失函数的特点：非负，且显然误分类点越少，或误分类点离超平面越近，损失函数值越小。直到没有误分类点时损失函数值为0。

感知机的学习策略就是在假设空间寻找 $w,b$ 让对该数据集的所有点，损失函数最小。

学习算法

优化损失函数的方法有很多，最简单是梯度下降。对感知机的损失函数，两个参数的梯度
$\nabla_w \ L(w,b) = \ -\sum \limits_{x_i \in M} y_i*x_i\\\nabla_b \ L(w,b) = \ -\sum \limits_{x_i \in M} y_i$

所以每次更新两个参数的迭代式子是 $w =w+h*y_i*x_i$ 和 $b= b+h*y_i$ ，其中 $h$ 是步长，范围在0到1之间，又称为学习率。

感知机学习算法的原始形式步骤：

选取初值 $w_0,b$
在训练集中选取数据 $(x_i,y_i)$
如果该数据是误分类点，即 $y_i(wx_i+b)\le 0$ ，更新参数
$w =w+h*y_i*x_i\\b= b+h*y_i$
转向2. ，直到没有误分类点，或者达到训练轮次。

上述步骤被称为感知机学习算法的原始形式。

感知机算法的对偶形式

和原始形式并无太大差别，但是可以加速训练。

对偶形式的思想：将 $w$ 和 $b$ 看作是实例和标签的线性组合。对于每一个样本 $(x_i,y_i)$ ，在更新过程中使用了 $n_i$ 次，即总次循环中，有 $n_i$ 次中将该样本作为了误分类点，故用它去更新参数。而一共有N个样本。

原始形式 $w =w+h*y_i*x_i$ 和 $b= b+h*y_i$ 就可以写成
$w=\sum\limits_{i=1} \limits^{N} n_i hy_ix_i \\b = \sum\limits_{i=1} \limits^{N} n_ihy_i$
则，感知机模型化为
$f(x)=sign(w\cdot x+b)=sign(\sum\limits_{i=1} \limits^{N} n_i hy_ix_i \cdot x + \sum\limits_{i=1} \limits^{N} n_ihy_i)$
学习目标变成了 $n_i$ 。

训练过程如下：

学习参数 $n=(n_1,n_2,....n_N)$ ，初始赋值全部是0
在数据集中选取数据 $(x_j,y_j)$
判断是不是误分类点，即 $y_j (\sum\limits_{i=1} \limits^{N} n_i hy_ix_i \cdot x_j + \sum\limits_{i=1} \limits^{N} n_ihy_j) \le 0$ ，如果是，更新 $n_i = n_i+1$
转至2，直到没有误分类点

可以从对偶形式的计算式子中看到，样本之间的计算都是 $x_i \cdot x_j$ ，其余计算都是N维向量的矩阵。其中N是样本个数，因此对偶形式适用于样本个数比特征空间的维数小的情况。

样本之间的內积计算，可以在一开始就计算存储为Gram矩阵，即 $G=[x_i \cdot x_j]_{N\times N}$ ，进行参数更新时之间查表即可，可以加快训练速度。

程序如下，使用的是CSV版本的mnist数据集。train.CSV是60000行，每一行的第一列是该图片的类别，第二列到第785列分别是该图片从上到下从左到右的位置处的像素值。
这个数据集网上到处都有，随便下载就行。

'''
@Author: muyi
@Date: 2020-01-04 21:02:27
@LastEditors  : muyi
@LastEditTime : 2020-01-05 23:22:41
@Description: 效果如下
在10000个样本中测试,4861个预测错误,正确率是0.513900
用时: 27.678825616836548
'''
import numpy as np
import time

def load_data(filename):
    '''
    加载mnist数据集，CSV格式，返回的是列表形式
    '''
    data=[]; label = []
    with open(filename) as f:
        for line in f.readlines():

            line = line.strip().split(",")  # csv用逗号分隔开
            if int(line[0])>=5:
                label.append(1)
            else:
                label.append(-1)
            
            data.append(list(map(int,line[1:]))) # 一行有785个，从第二个到最后一个是图像的像素值 在0-255之间
    return data,label


def perceptron(data,labels,epoches= 50,h=0.001):
    '''
    输入数据和标签进行训练，迭代次数epoch和步长h
    '''
    data = np.array(data)  # 也可以进行归一化，除以255，随便
    labels = np.array(labels)
    m, n = data.shape  # 数据集的行列数 m个数据 每个数据n列，即n个特征 mnist中是784
    
    #w = np.zeros(shape=(1,n))
    mu,sigma=0,0.1 #均值与标准差 w初值随机生成
    w = np.random.normal(mu,sigma,n).reshape(1,n)
    b = 0

    for epoch in range(epoches):  # 每一个批次 全部样本送入训练
        for i in range(m): # 遍历m个样本
            xi = data[i]
            yi = labels[i]
            if yi*(np.dot(w,xi)[0]+b) <=0 : # 判断是不是误分类点
                w = w + h*yi*xi
                b = b + h*yi
    print("training end")
    return w,b

def calculate_acc(data,label, w,b):  
    '''
    计算测试集上的准确率
    输入测试集和训练的参数w,b
    '''
    data = np.array(data)
    labels = np.array(label)
    m, n = data.shape
    count = 0
    for i in range(m):
        xi = data[i]
        yi = label[i]
        prediction = yi*(np.dot(xi,b)[0] + b)
        if prediction <0:
            count += 1  # 错误的样本数目
    acc = 1- count / m
    print("在%d个样本中测试,%d个预测错误,正确率是%f"%(m,count,acc))
    return acc

def dual_perce(data,label,epoches=10,h=0.001):
    data = np.array(data)/255.
    label = np.array(label)
    m, n = data.shape  # 数据集的行列数 m个数据 每个数据n列，即n个特征 mnist中是784
    G = np.zeros(shape=(m,m),dtype=np.uint8)  # Gram矩阵是样本的互相內积
    # 计算Gram矩阵
    for i in range(m):
        for j in range(m):
            tmp = np.dot(data[i].reshape(1,784),data[j])[0]
            G[i][j] = tmp
            G[j][i] = tmp

    # 初始化参数a,其中a[i]表示第i个样本的更新次数
    a = np.zeros(shape=(m,1))
    
    for epoch in range(epoches):
        for j in range(m): # 遍历每一个样本
            xj = data[j]
            yj = label[j]
            # 判断是不是误分类点
            w = sum([a[i]*h*label[i]*G[i][j] for i in range(m)  ])
            b = sum([a[i]*h*yj for i in range(m)])
            if yj*(w+b) <=0 :
                a[j] = a[j]+1
    w = sum([a[i]*h*label[i]*data[i] for i in range(m)])
    b = sum([a[i]*h*label[i] for i in range(m)])
    return w,b

if __name__ == "__main__":
    strat = time.time()
    train_x,train_y= load_data(r"mnist_train.csv")
    w, b = perceptron(train_x,train_y,50,h=0.0001)
    # w, b = dual_perce(train_x,train_y,10,h=0.0001) 对偶形式，事实上这个会特别慢，
    #因为mnist特征空间是784维度，而样本数60000，对偶形式在样本数小于特征空间维数时才有效果 

    test_x,test_y= load_data(r"mnist_test.csv")
    calculate_acc(test_x,test_y,w,b)
    end = time.time()
    print("用时:",end-strat)

欢迎勘误。