统计学习方法---感知机

《统计学习方法》系列笔记的第一篇，对应原著第二章。大量引用原著讲解，加入了自己的理解。对书中算法采用Python实现，并用Matplotlib可视化了动画出来。

概念

感知机是二分类模型，输入实例的特征向量，输出实例的±类别。

感知机模型

定义

假设输入空间是，输出空间是，x和y分属这两个空间，那么由输入空间到输出空间的如下函数：

称为感知机。其中，w和b称为感知机模型参数，叫做权值或权值向量，叫做偏置，w·x表示向量w和x的内积。sign是一个函数：

感知机的几何解释是，线性方程

将特征空间划分为正负两个部分：

这个平面（2维时退化为直线）称为分离超平面。

感知机学习策略

数据集的线性可分性

定义

给定数据集

其中如果存在某个超平面S

能够完全正确地将正负实例点全部分割开来，则称T线性可分，否则称T线性不可分。

感知机学习策略

假定数据集线性可分，我们希望找到一个合理的损失函数。

一个朴素的想法是采用误分类点的总数，但是这样的损失函数不是参数w，b的连续可导函数，不可导自然不能把握函数的变化，也就不易优化（不知道什么时候该终止训练，或终止的时机不是最优的）。

另一个想法是选择所有误分类点到超平面S的总距离。为此，先定义点x0到平面S的距离：

分母是w的L2范数，所谓L2范数，指的是向量各元素的平方和然后求平方根（长度）。这个式子很好理解，回忆中学学过的点到平面的距离：

此处的点到超平面S的距离的几何意义就是上述距离在多维空间的推广。

又因为，如果点i被误分类，一定有

成立，所以我们去掉了绝对值符号，得到误分类点到超平面S的距离公式：

假设所有误分类点构成集合M，那么所有误分类点到超平面S的总距离为

分母作用不大，反正一定是正的，不考虑分母，就得到了感知机学习的损失函数：

感知机学习算法

原始形式

感知机学习算法是对以下最优化问题的算法：

感知机学习算法是误分类驱动的，先随机选取一个超平面，然后用梯度下降法不断极小化上述损失函数。损失函数的梯度由：

给出。所谓梯度，是一个向量，指向的是标量场增长最快的方向，长度是最大变化率。所谓标量场，指的是空间中任意一个点的属性都可以用一个标量表示的场（个人理解该标量为函数的输出）。

随机选一个误分类点i，对参数w，b进行更新：

上式是学习率。损失函数的参数加上梯度上升的反方向，于是就梯度下降了。所以，上述迭代可以使损失函数不断减小，直到为0。于是得到了原始形式的感知机学习算法：

对于此算法，使用下面的例子作为测试数据： 

       

         

                 

      

给出Python实现和可视化代码如下：

感知机算法代码

终于到了最激动人心的时刻了，有了上述知识，就可以完美地可视化这个简单的算法：

1. # -*- coding:utf-8 -*-
2. 
   
3. import copy
4. from matplotlib import pyplot as plt
5. from matplotlib import animation
6.  
7. training_set = [[(3, 3), 1], [(4, 3), 1], [(1, 1), -1]]
8. w = [0, 0]
9. b = 0
10. history = []
11.  
12.  
13. def update(item):
14.     """
15.     update parameters using stochastic gradient descent
16.     :param item: an item which is classified into wrong class
17.     :return: nothing
18.     """
19.     global w, b, history
20.     w[0] += 1 * item[1] * item[0][0]
21.     w[1] += 1 * item[1] * item[0][1]
22.     b += 1 * item[1]
23.     print w, b
24.     history.append([copy.copy(w), b])
25.     # you can uncomment this line to check the process of stochastic gradient descent
26.  
27.  
28. def cal(item):
29.     """
30.     calculate the functional distance between 'item' an the dicision surface. output yi(w*xi+b).
31.     :param item:
32.     :return:
33.     """
34.     res = 0
35.     for i in range(len(item[0])):
36.         res += item[0][i] * w[i]
37.     res += b
38.     res *= item[1]
39.     return res
40.  
41.  
42. def check():
43.     """
44.     check if the hyperplane can classify the examples correctly
45.     :return: true if it can
46.     """
47.     flag = False
48.     for item in training_set:
49.         if cal(item) <= 0:
50.             flag = True
51.             update(item)
52.     # draw a graph to show the process
53.     if not flag:
54.         print "RESULT: w: " + str(w) + " b: " + str(b)
55.     return flag
56.  
57.  
58. if __name__ == "__main__":
59.     for i in range(1000):
60.         if not check(): break
61.  
62.     # first set up the figure, the axis, and the plot element we want to animate
63.     fig = plt.figure()
64.     ax = plt.axes(xlim=(0, 2), ylim=(-2, 2))
65.     line, = ax.plot([], [], 'g', lw=2)
66.     label = ax.text([], [], '')
67.  
68.     # initialization function: plot the background of each frame
69.     def init():
70.         line.set_data([], [])
71.         x, y, x_, y_ = [], [], [], []
72.         for p in training_set:
73.             if p[1] > 0:
74.                 x.append(p[0][0])
75.                 y.append(p[0][1])
76.             else:
77.                 x_.append(p[0][0])
78.                 y_.append(p[0][1])
79.  
80.         plt.plot(x, y, 'bo', x_, y_, 'rx')
81.         plt.axis([-6, 6, -6, 6])
82.         plt.grid(True)
83.         plt.xlabel('x')
84.         plt.ylabel('y')
85.         plt.title('Perceptron Algorithm ()')
86.         return line, label
87.  
88.  
89.     # animation function.  this is called sequentially
90.     def animate(i):
91.         global history, ax, line, label
92.  
93.         w = history[i][0]
94.         b = history[i][1]
95.         if w[1] == 0: return line, label
96.         x1 = -7
97.         y1 = -(b + w[0] * x1) / w[1]
98.         x2 = 7
99.         y2 = -(b + w[0] * x2) / w[1]
100.         line.set_data([x1, x2], [y1, y2])
101.         x1 = 0
102.         y1 = -(b + w[0] * x1) / w[1]
103.         label.set_text(history[i])
104.         label.set_position([x1, y1])
105.         return line, label
106.  
107.     # call the animator.  blit=true means only re-draw the parts that have changed.
108.     print history
109.     anim = animation.FuncAnimation(fig, animate, init_func=init, frames=len(history), interval=1000, repeat=True,
110.                                    blit=True)
111.     plt.show()
112.     anim.save('perceptron.gif', fps=2, writer='imagemagick')

可视化

可见超平面被误分类点所吸引，朝着它移动，使得两者距离逐步减小，直到正确分类为止。通过这个动画，是不是对感知机的梯度下降算法有了更直观的感悟呢？

算法的收敛性

记输入向量加进常数1的拓充形式，其最大长度为，记感知机的参数向量，设满足条件的超平面可以将数据集完全正确地分类，定义最小值伽马：

则误分类次数k满足：

证明请参考《统计学习方法》P31。

感知机学习算法的对偶形式

对偶指的是，将w和b表示为测试数据i的线性组合形式，通过求解系数得到w和b。具体说来，如果对误分类点i逐步修改wb修改了n次，则w，b关于i的增量分别为，这里，则最终求解到的参数分别表示为：

于是有算法2.2：

感知机对偶算法代码

涉及到比较多的矩阵计算，于是用NumPy比较多：

1. # -*- coding:utf-8 -*-
2. 
   
3. import numpy as np
4. from matplotlib import pyplot as plt
5. from matplotlib import animation
6.  
7. # An example in that book, the training set and parameters' sizes are fixed
8. training_set = np.array([[[3, 3], 1], [[4, 3], 1], [[1, 1], -1]])
9.  
10. a = np.zeros(len(training_set), np.float)
11. b = 0.0
12. Gram = None
13. y = np.array(training_set[:, 1])
14. x = np.empty((len(training_set), 2), np.float)
15. for i in range(len(training_set)):
16.     x[i] = training_set[i][0]
17. history = []
18.  
19. def cal_gram():
20.     """
21.     calculate the Gram matrix
22.     :return:
23.     """
24.     g = np.empty((len(training_set), len(training_set)), np.int)
25.     for i in range(len(training_set)):
26.         for j in range(len(training_set)):
27.             g[i][j] = np.dot(training_set[i][0], training_set[j][0])
28.     return g
29.  
30.  
31. def update(i):
32.     """
33.     update parameters using stochastic gradient descent
34.     :param i:
35.     :return:
36.     """
37.     global a, b
38.     a[i] += 1
39.     b = b + y[i]
40.     history.append([np.dot(a * y, x), b])
41.     # print a, b # you can uncomment this line to check the process of stochastic gradient descent
42.  
43.  
44. # calculate the judge condition
45. def cal(i):
46.     global a, b, x, y
47.  
48.     res = np.dot(a * y, Gram[i])
49.     res = (res + b) * y[i]
50.     return res
51.  
52.  
53. # check if the hyperplane can classify the examples correctly
54. def check():
55.     global a, b, x, y
56.     flag = False
57.     for i in range(len(training_set)):
58.         if cal(i) <= 0:
59.             flag = True
60.             update(i)
61.     if not flag:
62.  
63.         w = np.dot(a * y, x)
64.         print "RESULT: w: " + str(w) + " b: " + str(b)
65.         return False
66.     return True
67.  
68.  
69. if __name__ == "__main__":
70.     Gram = cal_gram()  # initialize the Gram matrix
71.     for i in range(1000):
72.         if not check(): break
73.  
74.     # draw an animation to show how it works, the data comes from history
75.     # first set up the figure, the axis, and the plot element we want to animate
76.     fig = plt.figure()
77.     ax = plt.axes(xlim=(0, 2), ylim=(-2, 2))
78.     line, = ax.plot([], [], 'g', lw=2)
79.     label = ax.text([], [], '')
80.  
81.     # initialization function: plot the background of each frame
82.     def init():
83.         line.set_data([], [])
84.         x, y, x_, y_ = [], [], [], []
85.         for p in training_set:
86.             if p[1] > 0:
87.                 x.append(p[0][0])
88.                 y.append(p[0][1])
89.             else:
90.                 x_.append(p[0][0])
91.                 y_.append(p[0][1])
92.  
93.         plt.plot(x, y, 'bo', x_, y_, 'rx')
94.         plt.axis([-6, 6, -6, 6])
95.         plt.grid(True)
96.         plt.xlabel('x')
97.         plt.ylabel('y')
98.         plt.title('Perceptron Algorithm 2 (www.hankcs.com)')
99.         return line, label
100.  
101.  
102.     # animation function.  this is called sequentially
103.     def animate(i):
104.         global history, ax, line, label
105.  
106.         w = history[i][0]
107.         b = history[i][1]
108.         if w[1] == 0: return line, label
109.         x1 = -7.0
110.         y1 = -(b + w[0] * x1) / w[1]
111.         x2 = 7.0
112.         y2 = -(b + w[0] * x2) / w[1]
113.         line.set_data([x1, x2], [y1, y2])
114.         x1 = 0.0
115.         y1 = -(b + w[0] * x1) / w[1]
116.         label.set_text(str(history[i][0]) + ' ' + str(b))
117.         label.set_position([x1, y1])
118.         return line, label
119.  
120.     # call the animator.  blit=true means only re-draw the parts that have changed.
121.     anim = animation.FuncAnimation(fig, animate, init_func=init, frames=len(history), interval=1000, repeat=True,
122.                                    blit=True)
123.     plt.show()
124.     # anim.save('perceptron2.gif', fps=2, writer='imagemagick')

可视化

与算法1的结果相同，我们也可以将数据集改一下：

training_set = np.array([[[3, 3], 1], [[4, 3], 1], [[1, 1], -1], [[5, 2], -1]])

会得到一个复杂一些的结果：

读后感

通过最简单的模型，学习到ML中的常用概念和常见流程。

另外本文只是个人笔记，服务于个人备忘用，对质量和后续不做保证。还是那句话，博客只做补充，要入门，还是得看经典著作。

Reference

文中部分代码参考了OldPanda的实现。