线性判别函数,感知器和松弛算法

线性判别函数是SVM、神经网络等技术的基础。本文主要内容是由线性判别函数得到的线性分类器、判别函数学习算法和感知机。线性判别函数要求问题是线性可分的,对于线性不可分的问题,可以将特征映射到高维空间再寻找分界面,再使用线性判别函数的广义形式进行分类。线性判别函数的学习过程可以转化为最优化损失函数的过程,即定义损失函数再使用梯度下降法进行学习。感知器会作为学习线性分类器的例子在本文进行讨论,其为多层神经网络的基础。关于感知器的内容主要是其定义及其学习算法。最后会讨论判别函数学习中的松弛算法。其中,感知器和松弛算法会在Python语言下进行实现。

线性判别函数

对于分类问题,在已有一定量样本的情况下(已知每个样本的特征及类别),若希望根据样本集学习一个分类器,使得对于一个未知类别的样本$x$能根据其的特征进行分类,一种做法是假设数据总体服从某个概率分布,即有概率分布$P(w_i|x)$(其中$i=1,2,…,c$, $w_i$ 表示第$i$个类别)存在,再从样本集中学习这些概率分布再进行分类(如参数估计和非参估计)。另一种方法是跳过这些概率推理的过程,直接假设最后得到的分类器满足一定的函数形式。线性分类器就属于后一种方法,其假设分类器的函数形式是线性判别函数(linear discriminant functions)。线性判别函数是输入$x$的线性(一次)函数,其表达式为:
$$
g(x) = w^tx + w_o
$$
其中,$w$是权重向量,$w_0$是偏差值。下面基于两类分类问题对线性判别函数进行讨论。在只有两类的样本集中,线性判别函数将样本分为两个类别,如果 $g(x)>0 $ ,则$x \in w_1$ ,否则$x \in w_2$。因$g(x) = 0$是$d$维(样本$x$的特征数)空间中的一个超平面,$g(x)$将$d$维空间分隔为两个空间,每个空间对应一个类别。对于位于分隔超平面上的任意点$x_1, x_2$,由于$g(x)$,故有$w^t(x_1 - x_2) = 0$,由于$x_1-x_2$是位于分隔超平面上向量,故$w^t$与分隔超平面上的任意向量正交,即$w_t$与分隔超平面垂直。对于空间中的任一点$x$,其可表示为$x = x_p + r \frac{w}{||w||}$,其中$x_p$是$x$在分隔超平面的投影点坐标,$r$为$x_p$到$x$距离,故有:
$$
g(x) = g(x_p+r \frac{w}{||w||}) = w_t(x_p + r\frac{w}{||w||}) + w_0 = r||w||
$$
对于$w \neq 0$ 有$r = g(x)/||w||$,故对于固定的$w$,$g(x)$与$x$点到分隔超平面的距离成正比

广义线性判别函数

线性判别函数使用超平面将特征空间分隔为两个部分,学习一个线性判别函数需要要求不同类别样本在特征空间下能被超平面区分(线性可分)。当样本不能满足线性可分时,不能求出上文提到的分隔超平面。然而,我们可以对线性判别函数形式进行推广,从而将样本特征映射到高维空间中,再在高维空间中寻找分隔超平面。对线性判别函数进行推广可得广义线性判别函数(Generalized Linear Discriminant Functions),其表达式为:
$$
g(x) = w^ty(x) + w_o
$$
其中,$y(x)$的作用是将$d$维向量$x$映射到$k$维空间,$w^t$权重向量的维数是$k$。例如,对于二维平面上的样本$x = (x_1, x_2)^t$,映射函数$y(x) = (x_1, x_2, x_1x_2)^t$可以将二维点映射到三维空间中。寻找合适的映射函数通常是一个困难的过程,SVM给出一种便捷的解决方式。

线性判别函数的学习算法:梯度下降法

知道了线性判别函数的形式后,那么下个问题便是如何获得线性判别函数,即如何进行学习。从线性判别函数的表达式来看,学习的目标既找到$w^t, w_0$的数值,使得其能区分样本集的类别。进行线性判别函数学习的一种方式是通过梯度下降法。梯度下降法的思想是将所有可能的$w^t,w_0 $值看成候选值,并对候选值进行打分(定义一个打分函数),假如能找到打分最高的候选值就意味着找到了最佳的$w^t, w_0$值。由于打分函数通常是候选值的损失函数,故找最佳权重变成最小化损失函数。找最佳候选值的过程使用了打分函数的梯度值,减少损失函数值意味则沿着损失函数的梯度下降,故称梯度下降法。

下面来形式化以上语言。首先,合并权重$w$和偏移量$w_0$使得$a = (w^t, w_0)^t$,再对样本进行预处理使得所有样本属于同一类别,操作为$y = y_c(y_1(x), y_2(x), …., y_k(x), 1)$,其中$y_c$表示样本$x$的类别,若属于第一类,则$y_c = 1$,否则为$y_c = -1$,这样求解广义线性判别方程就求$a$使得$g(x) = a^ty \geq 0$对所有样本成立。

损失函数是特征向量权重值的函数,记为$J(a)$。对损失函数求导得到梯度 $ \nabla J(a)$,一般地对于$a_0$,沿其梯度 $ \nabla J(a_0)$反方向进行移动能减少损失函数值$J(a_0)$。故比$a_0$更好的选择是向$\nabla J(a_0)$方向移动,将移动的过程乘上系数 $ \eta $,再不断迭代即是梯度下降法,迭代公式为:
$$
a(k+1) = a(k) - \eta \nabla J(a(k))
$$
其中,$a(k)$表示为第$k$迭代的权重值。梯度下降法的伪代码如下:

1
2
3
4
5
6
# compute gradient descent
Initialize a, k = 0, max_epochs, gradient function G(a, data), threshold \theta
while G(a) >= \theta or k < max_epochs:
a_{k+1} = a_k - \eta * G(a_k, data)
k += 1
return a_k

在进行迭代学习的过程中,梯度更新的过程可以在样本集上进行-批量学习(Batch Training),也可以对单个样本进行随机学习(Stochastic Training)。

感知器

感知器是Frank Rosenblatt在1957年发明的一种人工神经网络,其可视为单层神经网络。感知器使用了线性判别函数进行样本判别,故其要求样本在特征空间中线性可分。感知器定义损失函数为:$J_p(a) = \sum (-a^t y’)$,其中$y’$表示所有被判别函数分类错误的样本。其梯度$\nabla J_p (a) = \sum (- y’)$。因此感知器的迭代方程为:$ a(k+1) = a(k) + \eta (y’) $,即迭代的过程是每次从全体样本集中找到错误分类的样本特征,再添加到权重值中。实现的时候可以每次迭代时可以只针对一个样本(随机梯度下降),即如遇到错误分类则更新:$a(k+1) = a(k) + \eta y’$。Novikoff(1962)证明如果训练集是线性分隔的,则此种每次更新一个样本的算法(随机学习)会在有限次迭代后收敛。

mlxtend 的源码给出很好的Python实现范例(见其说明文档12 classifier.Perceptron章节),按其框架在下面进行实现。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
import numpy as np
# Rosenblatt Perceptron
# classifier = Perceptron(eta = step_len, epochs = max_iteration)
# classifier.train(X, y)
# make sure X = [nsample X nfeature], y = [nsample X 1]
# y only have two class label {1, -1}
class Perceptron(object):
def __init__(self, eta=0.01, epochs=50):
self.eta = eta # eta should > 0
self.epochs = epochs
# w_ = [(nfeature + 1) X 1]
# 随机梯度下降
def increment_train(self, X, y, w_ = None, eta=None, epochs=None):
# init weight
self.w_ = np.zeros(1 + X.shape[1]) if w_ is None \
else np.array(w_)
self.eta = self.eta if eta is None else eta
self.epochs = self.epochs if epochs is None else epochs
for _ in range(self.epochs):
#update weights
error = 0
for xi, target in zip(X, y):
pre = self.predict(xi)
if target != pre:
error += 1
self.w_[0] += target * self.eta
self.w_[1:] += target * self.eta * xi
print _, error
if error == 0:
break
# 批量梯度下降
def batch_train(self, X, y, w_ = None, eta=None, epochs=None):
# init weight
self.w_ = np.zeros(1 + X.shape[1]) if w_ is None \
else np.array(w_)
self.eta = self.eta if eta is None else eta
self.epochs = self.epochs if epochs is None else epochs
for _ in range(self.epochs):
#update weights
error = 0
update_0, update_ = 0, 0
for xi, target in zip(X, y):
pre = self.predict(xi)
if target != pre:
error += 1
update_0 += target
update_ += target * xi
print _, error
if error == 0:
break
else:
self.w_[0] += self.eta * update_0
self.w_[1:] += self.eta * update_
# net ouput, w^t * X
def net_activation(self, X):
return np.dot(X, self.w_[1:]) + self.w_[0]
# transfer or predict 1 if net > 0
def predict(self, X):
return np.where(self.net_activation(X) >= 0.0, 1, -1)
def get_predict_rate(self, X, y):
return sum(self.predict(X) == y) * 1.0 / X.shape[0]

下面使用使用Iris 数据集进行测试。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import pandas as pd
# read data
df = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', header=None)
y = df.iloc[0:100, 4].values
y = np.where(y == 'Iris-setosa', -1, 1)
X = df.iloc[0:100, [0,2]].values
%matplotlib inline
import matplotlib.pyplot as plt
from mlxtend.plotting import plot_decision_regions
ppn = Perceptron(epochs=10, eta=0.1)
ppn.increment_train(X, y)
#ppn.batch_train(X, y)
print('Weights: %s, predict rate:%s' % (ppn.w_, ppn.get_predict_rate(X, y)))
plot_decision_regions(X, y, clf=ppn)
plt.title('Perceptron')
plt.xlabel('sepal length [cm]')
plt.ylabel('petal length [cm]')
plt.show()

在Iris 数据集下使用随机梯度下降的权重是 $w = [-0.2 , -0.34 ,0.91]^t$, 使用批量梯度下降得到的权重是 $w = [ -5.3 ,-12.81,33.18]^t$,可见后者的数值远大于前者(批量梯度下降的过程中每次迭代的权重值是所有所有错分样本特征的和),但是两者分类图像几乎相当,可见权重数乘的结果不影响分类结果。

prceptron_plot

上图只迭代了10次就得到线性分界面。感知器的算法十分简洁,但实际应用会存在许多问题:(1)当数据线性可分时感知器得到的分类界面不唯一确定,其与初始权重有关;(2)收敛速度可能会很慢,但数据可选分类界面十分狭窄时迭代次数会很大;(3)若数据线性不可分,算法不收敛并分类界面会循环移动,循环的周期可能很大,故可能难以察觉。

感知器的学习过程需要设置合适的步长$\eta$ 值,当样本的每个特征的数值差异较大时设置合适的步长是困难的。一种改善该问题的预处理方法是对特征分量进行标准,如将每个分量映射到$[0,1]$区间,或将其转换为$(0, 1)$正态分布。

1
2
3
4
5
6
7
8
9
# normalize
X_std = np.copy(X)
X_std[:,0] = (X[:,0] - X[:,0].mean()) / X[:,0].std()
X_std[:,1] = (X[:,1] - X[:,1].mean()) / X[:,1].std()
# scaling
X_scale = np.cory(X)
X_scale[:,0] = (X[:,0] - X[:,0].min()) / (X[:,0].max() - X[:,0].min())
X_scale[:,1] = (X[:,1] - X[:,1].min()) / (X[:,1].max() - X[:,1].min())

松弛学习算法

感知器的损失函数是$J(a) = \sum -a^ty’$,由于$g(y’) = a^ty’= r ||a||$,其中$r$是该错分点$y’$到分解面的距离。对于固定的$y’$,将该函数看为权重值的函数,则有$g(a) = a^ty’= r_2||y’||$,其中$r_2$为在权重空间下权重点 $a$ 到分解面的距离。故有 $r_2 = g(a) / ||y’|| = a^ty’ / ||y’||$。考虑将此值所谓损失度量,则有$J(a) = r_2^2 = \frac{(a^ty’)^2}{ 2||y’||^2}$ 。不失一般性将判别函数$a^ty > 0$ 转换为 $a^ty > b$ ,其中$b\geq 0$,表示当线性组合值大一定值时才判为正类,这样权重到平面的距离即为 $(a^ty’ - b)/ ||y’||$,损失函数即为$J(a) = \sum (r_2)^2 = \sum \frac{(a^ty’ - b)^2}{ 2||y’||^2}$,求梯度得到$ \nabla J(a) = \sum \frac{(a^ty’ - b)y’}{ ||y’||^2}$ ,该梯度即为$a$到分界面$a^ty = b$的距离!这样权重迭代更新规则变成$a(k+1)= a(k) - \eta \nabla J(a{k}) = a(k) - \eta \sum \frac{(a(k)^ty’ - b)y’}{ ||y’||^2}$,这就是松弛算法。因为每次对权重的更新都使得将$y’$(对随机梯度下降来说的单个样本)引起的错分进行”松弛“。

下面是Python代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
class Relaxation(object):
def __init__(self, eta=0.01, epochs=50, b=0.1):
self.eta = eta # eta should > 0
self.epochs = epochs
self.b = b
def train(self, X, y, w_ = None):
# init weight
self.w_ = np.zeros(1 + X.shape[1]) if w_ is None \
else np.array(w_)
for _ in range(self.epochs):
#update weights
error = 0
for xi, target in zip(X, y):
gx = np.dot(xi, self.w_[1:]) + self.w_[0]
pre = 1 if gx > self.b else -1
# print pre, target
if target != pre:
error += 1
update = self.eta * (gx - self.b) / (sum(xi**2) + 1)
self.w_[0] -= update * target
self.w_[1:] -= update * xi * target
print _, error, self.w_
if error == 0:
break
def predict(self, X):
return np.where((np.dot(X, self.w_[1:]) + self.w_[0]) > self.b, 1, -1)
def predict_rate(self, X, y):
return sum(self.predict(X) == y) * 1.0 / X.shape[0]
X_std = np.copy(X)
X_std[:,0] = (X[:,0] - X[:,0].mean()) / X[:,0].std()
X_std[:,1] = (X[:,1] - X[:,1].mean()) / X[:,1].std()
rlx = Relaxation(epochs=20, eta=0.1, b=0.1)
w_ = rlx.train(X, y)
print('Weights: %s, predict rate:%s' % (rlx.w_, rlx.predict_rate(X, y)))
plot_decision_regions(X, y, clf=rlx)
plt.title('Relaxation')
plt.xlabel('sepal length [cm]')
plt.ylabel('petal length [cm]')
plt.show()

实际使用的时候需要注意不能同时初始化$b=0,w=0$,此时权重值在迭代的过程中不会改变。小样本下测试可以看出松弛算法在实际应用的时候没有感知器效果好,即使在对样本进行标准化后在较少的迭代次数下没有得到线性分界面。

prceptron_plot

总结

通过学习线性判别函数得到的分类器可以很好地解决线性可分的分类问题。感知器是十分高效便捷的算法,其通过针对分类错误的样本的特征信息来更新权重。实际实现线性判别函数学习时需要考虑好步长,权重等值的初始化。还可以对样本进行标准化以得到更好的学习、分类效果。

参考

[1] Pattern classification: https://www.amazon.com/Pattern-Classification-Pt-1-Richard-Duda/dp/0471056693

[2] Perceptron_wikipedia: https://www.wikiwand.com/zh-hans/%E6%84%9F%E7%9F%A5%E5%99%A8

[3] mlxtend :https://sebastianraschka.com/pdf/software/mlxtend-latest.pdf

[4] The element of statistical learning: https://www.amazon.com/Elements-Statistical-Learning-Prediction-Statistics/dp/0387848576