机器学习(2) 感知机原理及实现

前言

        在上一篇博文机器学习(1)泛化误差上界的实现及分析中,分析了评价模型迁移学习能力的指标之一泛化误差。由Hoeffding不等式表明,在假设空间有限、样本空间有限的情况下,泛化误差是依据概率服从于泛化误差上界的,也就是说可以预估了解一个模型的可能犯错的能力。
        在这篇博文中,将给出第一个二分类模型——感知机模型,并介绍在点集线性可分的条件下,可以采用随机梯度下降算法以找到一个可以将正负点集分离开的超平面。最后,给出感知机模型随机梯度下降方法及其对偶方法的C++实现。

参考文献:《统计学习方法》,李航。

感知机模型

        感知机模型是最简单的、用于区分两类线性可分数据集的模型。数据集 X X X中所有的点都由每个特征的不同值唯一确定,每个数据的维数是特征的数量,记为 n n n,如果有 N N N个数据,则每个数据都记为 x i = ( x i ( 1 ) , x i ( 2 ) , x i ( 3 ) , . . . , x i ( n ) ) T , i ∈ [ 1 , N ] , x i ∈ X x_i=(x^{(1)}_i, x^{(2)}_i, x^{(3)}_i, ..., x^{(n)}_i)^{T},i\in[1,N],x_i\in{X} xi=(xi(1),xi(2),xi(3),...,xi(n))T,i[1,N],xiX;每个数据都有一个对应的类别标签,记为 y i y_i yi,我们将两种不同的数据集对应的类别分别定义为正类和负类,正类用 + 1 +1 +1来表示,负类用 − 1 -1 1表示,即 y i ∈ Y = { + 1 , − 1 } , i ∈ [ 1 , N ] y_i\in{Y}=\{+1, -1\},i\in[1,N] yiY={+1,1},i[1,N]。如果可以找到一个超平面,将正类的数据点全部划分至超平面的一侧, 而负类的数据点全部被划分至超平面的另一侧,我们称该数据集是线性可分的。
        为了更好理解数据集的概念,我们给出一个例子,希望用感知机模型来区分资深程序员和新手程序员。规定每个程序员都有两个特征:代码量、工资。这里有三个程序员, x 1 = { 3 , 3 } , x 2 = { 3 , 4 } , x 3 = { 1 , 1 } x_1=\{3,3\},x_2=\{3,4\},x_3=\{1,1\} x1={3,3},x2={3,4},x3={1,1},如下图所示。因为资深程序员代码写的多、工资挣得多,数据集中在二维平面的右上角;新手程序员代码写得少、工资挣得少,数据集中在二维平面的左下角,这种情况下就会存在一条直线将正类划分在直线上面,且将负类数据划分在直线下面。我们称这条直线为二维空间中的超平面(hyperplane)。

二维数据
        类似的,一维空间上的超平面是一个点,三维空间中的超平面就是传统意义上的平面,四维空间的超平面是一个三维组合,……,n维空间的超平面直观表现为n-1维的形象。不过大家大可不必费劲心思去想更多维度的超平面长成什么样子,因为机器学习中的每一个维度是一个特征,用于描述数据的特点,特点可以是一维的(代码量),二维的(代码量,工资),三维的(代码量,工资,发量),四维的(代码量,工资,发量,从业年龄),五维的(代码量,工资,发量,从业年龄,掌握算法数量),等等。给定维数之后,每个样本就是一个 n n n维特征组成的向量(向量指列向量),有N个数据就有N个向量。在 x = ( x ( 1 ) , x ( 2 ) , x ( 3 ) , . . . , x ( n ) ) T \bm x=(x^{(1)},x^{(2)}, x^{(3)},..., x^{(n)})^{T} x=(x(1),x(2),x(3),...,x(n))T上定义超平面为 w ⋅ x + b = 0 \bm w\cdot \bm x+b=0 wx+b=0,其中 w = ( w ( 1 ) , w ( 2 ) , w ( 3 ) , . . . , w ( n ) ) T \bm w=(w^{(1)},w^{(2)}, w^{(3)},..., w^{(n)})^{T} w=(w(1),w(2),w(3),...,w(n))T同样也是n维向量, w ⋅ x \bm w\cdot \bm x wx是两向量的内积。
        因此,感知机模型的求解结果就是找到一个超平面使得正负点集可以在超平面两侧分隔开。假设数据集是线性可分的,采用梯度下降法进行优求解;如果数据集不满足线性可分性,应采用其他学习策略。特别指出的是,如上图所示,正类和负类之间的空隙可能很大,因此只需要稍微调整一下斜率或者截距,就可以得到一个全新的符合要求的超平面,这也就是说,线性可分的数据集经过一定调整,一定可以在有限次数内找到一个合适的超平面使得点集在超平面两侧;同时,这样的超平面是不唯一的,会根据初始参数选择、权值调整的策略不同而得到不同的超平面。

感知机损失函数

        当感知机构成之后,就有了初始的超平面,这个超平面可能很好,让正负点集分离到平面的两侧;也可能是一个很差的超平面,让正集负集交织在了超平面的某一侧。为了表示超平面的分类性能好坏,引入损失函数的概念。损失函数满足这样的条件:

  1. 如果某个数据被正确分类,损失函数的数值为0。
  2. 如果某个数据被错误分类,损失函数的数值增加。并且分类出错的越离谱,损失函数的数值就越大;分类错误的数据越多,损失函数的数值也应该越大。

        因此,定义损失函数为所有分类错误点与超平面距离之和。点到线的距离是
d = 1 ∣ ∣ w ∣ ∣ ∣ w ⋅ x + b ∣ (1) d=\frac{1}{||\bm w||}|\bm w\cdot \bm x+b|\tag{1} d=w1wx+b(1)
        这里 ∣ ∣ w ∣ ∣ ||\bm w|| w是参数向量 w \bm w w L 2 L_2 L2范数。距离公式的证明参见我的另一篇博文数学基础知识(1) 点到超平面距离。为方便损失函数求导以实现梯度下降修改超平面的两个参数 w \bm w w b \bm b b,在损失函数中将距离中的正参数 1 ∣ ∣ w ∣ ∣ \frac{1}{||\bm w||} w1略去。
        为使输出标签为 { + 1 , − 1 } \{+1, -1\} {+1,1},选用sign函数获得用感知机输出的预测类别。记预测的类别为 p r e d _ y i pred\_y_i pred_yi有:
p r e d _ y i = s i g n ( f ( x ) ) = s i g n ( w ⋅ x + b ) = { + 1 , i f   w ⋅ x + b > 0 − 1 , i f   w ⋅ x + b ≤ 0 (2) pred\_y_i=sign(f(\bm x))=sign(\bm w\cdot \bm x+b)=\begin{dcases}&+1,if \space \bm w\cdot \bm x+b>0\\&-1,if\space \bm w\cdot\bm x +b\le0\end{dcases}\tag{2} pred_yi=sign(f(x))=sign(wx+b)={+1,if wx+b>01,if wx+b0(2)
        当预测成功时, p r e d _ y i = y i pred\_y_i=y_i pred_yi=yi,即 y i ( w ⋅ x + b ) > 0 y_i(\bm w\cdot \bm x+b)>0 yi(wx+b)>0
        而预测失败时,有 p r e d _ y i = − y i pred\_y_i=-y_i pred_yi=yi,即 y i ( w ⋅ x + b ) ≤ 0 y_i(\bm w\cdot \bm x+b)\leq0 yi(wx+b)0。由于我们需要找到一个超平面使得正集、负集全部处于超平面的两侧,当出现有数据点在超平面上的时候,依然认定当前的感知机平面没有收敛,需要继续进行迭代。
        因此,对于分类错误的数据点, ∣ w ⋅ x i + b ∣ = − y i ( w ⋅ x i + b ) |\bm w\cdot \bm {x_i}+b| = -y_i(\bm w\cdot \bm {x_i}+b) wxi+b=yi(wxi+b)。损失函数记为:
L ( w , b ) = − ∑ x i ∈ M y i ( w ⋅ x i + b ) (3) L(\bm w,\bm b)=-\sum\limits_{\bm {x_i}\in M}y_i(\bm w\cdot \bm {x_i}+b)\tag{3} L(w,b)=xiMyi(wxi+b)(3)
         M M M为所有误分类的点集合。

随机梯度下降法

        梯度下降法是机器学习过程中常用的一种优化损失函数的方法。想要了解梯度下降法,首先要感受一下什么是梯度。下面这张图像是一个峰值函数的图像,用Python绘制而成,是Matlab里面的peaks函数,我认为这张函数图像最适合于解释梯度下降法(而本文的误差函数是相反倒锥字形,选用这个图像仅仅是为了利用山峰的概念解释下山),所以在此展示出来。该函数的表达式为

  • z = 3 ( 1 − x ) 2 e − x 2 − ( y + 1 ) 2 − 10 ( x 5 − x 2 − y 5 ) e − x 2 − y 2 − 1 3 e − ( x + 1 ) 2 − y 2 \displaystyle{z=3(1-x)^2e^{-x^2-(y+1)^2}-10(\frac{x}{5}-x^2-y^5)e^{-x^2-y^2}-\frac{1}{3}e^{-(x+1)^2-y^2}} z=3(1x)2ex2(y+1)210(5xx2y5)ex2y231e(x+1)2y2

        图像的绘制方法参见我的另一篇博文Python数据分析(01) 绘制三维曲面图像。有了这张图像,我们就能够直观的感受到什么是梯度下降法了。
peaks

山峰函数全貌

山峰函数的正半轴部分

山峰函数的正半轴部分

        思来想去,我琢磨出了这样一句话:假设我们攀爬在一座由误差堆积而成的山峰之上,山峰之下是正确坦荡的平原。我们本应走在坦荡的平原大道之上,却误打误撞错攀上了误差的高峰。为了尽快使生活回归康庄大道,我们不得不选择一条此刻最快的下山之路。这座山峰是由一次次错误的选择堆砌而成,矫正过去犯下的错误,而又不至于矫枉过才是我们当下的选择。 我用了一段形而上的语言来总结随机梯度下降法的思路,下面我解释一下刚刚的这段话。

我们攀爬在一座由误差堆积而成的山峰之上,山峰之下是正确坦荡的平原。我们本应走在坦荡的平原大道之上,却误打误撞错攀上了误差的高峰。

        根据损失函数的定义,损失函数的数值是由逐个误差点积累而成的数值,我们将这种积累过程看做造山一样,认为误差会形成一座大大山;而当误差值为0的时候,就是感知机模型收敛的时候,产生了一个超平面使得正类样本和负类样本均匀的各自分布在了超平面的两侧,这就是所说的“平原”。

为了尽快使生活回归康庄大道,我们不得不选择一条此刻最快的下山之路。

        对于之前的损失函数 L ( w , b ) = − ∑ x i ∈ M y i ( w ⋅ x i + b ) L(\bm w,\bm b)=-\sum\limits_{\bm {x_i}\in M}y_i(\bm w\cdot \bm {x_i}+b) L(w,b)=xiMyi(wxi+b)而言,样本点已经给定的条件下,视 x i , y i \bm {x_i},y_i xi,yi为常量,造成损失函数增大减小的是损失函数中变量的 w , b \bm w,b w,b。这时用到了函数梯度的概念。关于梯度,参见我的另一篇博客数学基础知识(2) 梯度和方向向量。简而言之,梯度是一个与自变量维数相同的向量,梯度向量的方向指向了在函数上的某点增加速率最快的方向,梯度向量中各个维度的数值是函数对该自变量在该点处的偏导数,梯度的模表示函数在该点变化率的最大值。 为了沿着最陡峭的山路迅速下山,我们需要沿着梯度的反方向移动(反之,为了尽快上山,则需要沿着梯度的方向移动)。

  • 损失函数对自变量 w = ( w ( 0 ) , w ( 1 ) , . . . , w ( n ) ) T \bm w=(w^{(0)}, w^{(1)}, ..., w^{(n)})^T w=(w(0),w(1),...,w(n))T的梯度为:
    ∇ w L ( w , b ) = ( − ∑ x i ∈ M y i x i ( 0 ) , − ∑ x i ∈ M y i x i ( 1 ) , . . . , − ∑ x i ∈ M y i x i ( n ) ) T = − ∑ x i ∈ M y i ( x i ( 0 ) , x i ( 1 ) , . . . , x i ( n ) ) T = − ∑ x i ∈ M y i x i (3) \begin{aligned} \nabla_{\bm w}L(\bm w,b)=&(-\sum\limits_{\bm {x_i}\in M}y_ix_i^{(0)},-\sum\limits_{\bm {x_i}\in M}y_ix_i^{(1)}, ...,-\sum\limits_{\bm {x_i}\in M}y_ix_i^{(n)})^T\\ =&-\sum\limits_{\bm {x_i}\in M}y_i(x_i^{(0)},x_i^{(1)},...,x_i^{(n)})^T\\ =&-\sum\limits_{\bm {x_i}\in M}y_i\bm {x_i} \end{aligned}\tag{3} wL(w,b)===(xiMyixi(0),xiMyixi(1),...,xiMyixi(n))TxiMyi(xi(0),xi(1),...,xi(n))TxiMyixi(3)
  • 损失函数对自变量 b b b的梯度为:
    ∇ b L ( w , b ) = − ∑ x i ∈ M y i (4) \nabla_bL(\bm w,b) = -\sum\limits_{\bm {x_i}\in M}y_i\tag{4} bL(w,b)=xiMyi(4)

这座山峰是由一次次错误的选择堆砌而成,矫正过去犯下的错误,而又不至于矫枉过才是我们当下的选择。

        损失函数基于误分类点构成,除了沿梯度下降从全局降到全局最优解之外,也可以通过每个误差点进行考虑。想要减小损失函数的数值至0,考虑减少误分类点的个数,也就是逐个调整误差点至正确分类点的情况。从局部的情况考虑时,更有可能陷入局部最优解,通过设定每个误差点学习速率的情况加以限制。通过随机选择误差点进行全局修正的方法,我们称之为随机梯度下降法。随机选取一个误差点对 w , b w,b w,b进行修正,由上述公式推导可知,每次应当沿着梯度的反方向进行修正,修正的速率为 η ∈ ( 0 , 1 ] \eta\in(0,1] η(0,1],称为学习率。每次调整的过程如下:
w ← w − η ∇ w L ( w , b ) ∣ x i = w + η y i x i (5) \bm w\gets \bm w-\eta\nabla_{\bm w} L(\bm w,b)|_{\bm {x_i}}=\bm w+\eta y_i\bm {x_i}\tag{5} wwηwL(w,b)xi=w+ηyixi(5)
b ← b − η ∇ b L ( w , b ) ∣ x i = b + η y i (6) b\gets b-\eta\nabla_bL(\bm w,b)|_{x_i}=b+\eta y_i\tag{6} bbηbL(w,b)xi=b+ηyi(6)

真正的误差函数

更贴近损失函数的图像,下低上高。

        最后这张图更加贴合损失函数的真正模样。在线性可分的条件下,当全部数据集依据其正负样本的标签划分在超平面的两侧,则存在损失函数的极小值(同样也是最小值)0,但当超平面与数据集的距离增大时,损失函数的数值则会递增,趋向无穷大。因此,此前的山峰图像只是为了方便大家理解上山下山的问题,与感知机的损失函数不符。

感知机学习算法

        经过上面的讨论,我们给出感知机学习算法的原始形式:

  • 输入为原始数据 x i = ( x i ( 1 ) , x i ( 2 ) , . . . , x i ( n ) ) T , y i , η \bm {x_i}=(x_i^{(1)}, x_i^{(2)},...,x_i^{(n)})^T, y_i,\eta xi=(xi(1),xi(2),...,xi(n))T,yi,η,其中 i ∈ [ 1 , N ] , x i ∈ R n , y i ∈ { + 1 , − 1 } , η ∈ ( 0 , 1 ] i\in[1,N],\bm {x_i}\in R^n,y_i\in\{+1, -1\},\eta\in(0,1] i[1,N],xiRn,yi{+1,1},η(0,1]
  • 输出为超平面的参数 w = ( w ( 1 ) , w ( 2 ) , . . . , w ( n ) ) T , b \bm w=(w^{(1)}, w^{(2)}, ..., w^{(n)})^T,b w=(w(1),w(2),...,w(n))T,b,感知机模型 f ( x ) = s i g n ( w ⋅ x + b ) f(\bm x)=sign(\bm w\cdot \bm x+b) f(x)=sign(wx+b)
  1. 初始化参数 w 0 , b \bm {w_0}, b w0,b,一般初始化为0。
  2. 逐个训练数据 x i , y i \bm {x_i},y_i xi,yi
  3. 如果 y i ( w ⋅ x i + b ) ≤ 0 y_i(\bm w\cdot \bm {x_i}+b)\leq0 yi(wxi+b)0, w ← w + η   y i x i \bm w\gets \bm w+\eta\ y_i\bm {x_i} ww+η yixi b ← b + η y i b\gets b+\eta y_i bb+ηyi
  4. 重复步骤2,3,直至没有数据点被误分类。 ■ \blacksquare

        可以证明,在输入数据是线性可分的条件下,经过有限次迭代,一定可以将感知机损失函数收敛至零。依据 w , b \bm w,b w,b的初始化值不同、误差点选择的顺序不同、学习速率不同,最后生成的感知机超平面也未必相同,对同一组数据,有无限个超平面满足条件。感知机算法的收敛性参见定理Novikoff,此处不做更多展开。

感知机学习算法对偶形式

        感知机学习算法的原始形式中,每个数据都可能被误分类,在误分类的情况下会对 w , b w,b w,b进行调整,调整的结果仅与被误分类的次数、原始数据、学习速率有关。因此,学习后的感知机模型为:
w = w 0 + ∑ i = 1 N α i y i x i (7) \bm w=\bm {w_0}+\sum\limits_{i=1}^{N} \alpha_iy_i\bm {x_i}\tag{7} w=w0+i=1Nαiyixi(7)
b = b 0 + ∑ i = 1 N α i y i (8) b=b_0+\sum\limits_{i=1}^{N}\alpha_iy_i\tag{8} b=b0+i=1Nαiyi(8)
        其中, α i \alpha_i αi是第 i i i个数据被调整的次数 n i n_i ni乘以学习率 η \eta η,即 α i = n i η \alpha_i=n_i\eta αi=niη。因此,只需要在计算过程中记录每个数据被调整的次数,即可在调整结束时生成感知机的分类超平面。
        下面给出感知机算法对偶形式的表述。

  • 输入为原始数据 x i = ( x i ( 1 ) , x i ( 2 ) , . . . , x i ( n ) ) T , y i , η \bm {x_i}=(x_i^{(1)}, x_i^{(2)},...,x_i^{(n)})^T, y_i,\eta xi=(xi(1),xi(2),...,xi(n))T,yi,η,其中 i ∈ [ 1 , N ] , x i ∈ R n , y i ∈ { + 1 , − 1 } , η ∈ ( 0 , 1 ] i\in[1,N],\bm {x_i}\in R^n,y_i\in\{+1, -1\},\eta\in(0,1] i[1,N],xiRn,yi{+1,1},η(0,1]
  • 输出为超平面的参数 α = ( α 1 , α 2 , . . . , α N ) T , b \bm \alpha = (\alpha_1, \alpha_2, ..., \alpha_N)^T,b α=(α1,α2,...,αN)T,b,感知机模型 f ( x ) = s i g n ( ∑ i = 1 N α i y i x i ⋅ x + b ) f(\bm x)=sign(\sum\limits_{i=1}^N\alpha_iy_i\bm {x_i}\cdot \bm x+b) f(x)=sign(i=1Nαiyixix+b)
  1. 初始化参数 α , b \bm \alpha, b α,b,一般初始化为0。
  2. 逐个训练 x i , y i \bm {x_i}, y_i xi,yi
  3. 如果 y i ( ∑ j = 1 N α j y j x j ⋅ x i + b ) ≤ 0 y_i(\sum\limits_{j=1}^N\alpha_jy_j\bm {x_j}\cdot \bm {x_i}+b)\leq 0 yi(j=1Nαjyjxjxi+b)0 α i ← α i + η \alpha_i\gets\alpha_i+\eta αiαi+η b ← b + η y i b\gets b+\eta y_i bb+ηyi
  4. 重复步骤2,3,直至没有误差点出现。
            特别的,方便实现 x i ⋅ x j \bm {x_i}\cdot \bm {x_j} xixj的内积,在运算开始之前预计算所有数据两两之间的内积,记为Gram矩阵(Gram matrix)。 G = [ x i ⋅ x j ] N × N \bm{G}=[\bm {x_i}\cdot \bm {x_j}]_{N\times N} G=[xixj]N×N

代码实现

        下面的代码实现了感知机学习的原始算法和对偶算法两种方法,并通过控制输标志位,输出每次迭代的结果。以最初程序员的数据为例,使用原始算法,输出如下内容:

WrongNum = 1. Data[1] judge Wrong. w: 3.0 3.0 b: 1.0
$w·x + b: 3.0x^{(1)}+3.0x^{(2)}+1.0$
WrongNum = 2. Data[3] judge Wrong. w: 2.0 2.0 b: 0.0
$w·x + b: 2.0x^{(1)}+2.0x^{(2)}+0.0$
WrongNum = 3. Data[3] judge Wrong. w: 1.0 1.0 b: -1.0
$w·x + b: 1.0x^{(1)}+1.0x^{(2)}+-1.0$
WrongNum = 4. Data[3] judge Wrong. w: 0.0 0.0 b: -2.0
$w·x + b: 0.0x^{(1)}+0.0x^{(2)}+-2.0$
WrongNum = 5. Data[1] judge Wrong. w: 3.0 3.0 b: -1.0
$w·x + b: 3.0x^{(1)}+3.0x^{(2)}+-1.0$
WrongNum = 6. Data[3] judge Wrong. w: 2.0 2.0 b: -2.0
$w·x + b: 2.0x^{(1)}+2.0x^{(2)}+-2.0$
WrongNum = 7. Data[3] judge Wrong. w: 1.0 1.0 b: -3.0
$w·x + b: 1.0x^{(1)}+1.0x^{(2)}+-3.0$
Justify to accurate. w: 1.0 1.0 b: -3.0
w·x + b:
$1.0x^{(1)}+1.0x^{(2)}+-3.0$

        经过7次迭代,原始问题收敛至超平面:
w ⋅ x + b = x ( 1 ) + x ( 2 ) − 3 w\cdot x+b=x^{(1)}+x^{(2)}-3 wx+b=x(1)+x(2)3
        感知机收敛到:
f ( x ) = s i g n ( x ( 1 ) + x ( 2 ) − 3 ) f(x)=sign(x^{(1)}+x^{(2)}-3) f(x)=sign(x(1)+x(2)3)
        同样是程序员的例子,如果使用感知机学习的对偶算法,则会输出如下内容:

WrongNum = 1. Data[1] judge Wrong. alpha: 1 0 0 b: 1
WrongNum = 2. Data[3] judge Wrong. alpha: 1 0 1 b: 0
WrongNum = 3. Data[3] judge Wrong. alpha: 1 0 2 b: -1
WrongNum = 4. Data[3] judge Wrong. alpha: 1 0 3 b: -2
WrongNum = 5. Data[1] judge Wrong. alpha: 2 0 3 b: -1
WrongNum = 6. Data[3] judge Wrong. alpha: 2 0 4 b: -2
WrongNum = 7. Data[3] judge Wrong. alpha: 2 0 5 b: -3
Justify to accurate. alpha: 2 0 5 b: -3
w·x + b:
$1.0x^{(1)}+1.0x^{(2)}+-3.0$

        经过7次迭代, α \alpha α收敛到:
α = ( α 1 , α 2 , α 3 ) T = ( 2 , 0 , 5 ) T \alpha=(\alpha_1, \alpha_2, \alpha_3)^T=(2,0,5)^T α=(α1,α2,α3)T=(2,0,5)T
        原始问题收敛至超平面:
w ⋅ x + b = x ( 1 ) + x ( 2 ) − 3 w\cdot x+b=x^{(1)}+x^{(2)}-3 wx+b=x(1)+x(2)3
        感知机收敛到:
f ( x ) = s i g n ( x ( 1 ) + x ( 2 ) − 3 ) f(x)=sign(x^{(1)}+x^{(2)}-3) f(x)=sign(x(1)+x(2)3)

        另外,感知机算法适用于正类负类样本线性可分的情况,需要满足正实例点集构成的凸壳与负实例点构成的凸壳不存在交集。这需要用到两组向量组线性无关的判定方法,实现方法和证明方法较为繁琐。在我的程序中,简单的设置一个迭代上限,如果在迭代周期内没有实现收敛,就简单地判定这组样本是线性不可分的,可以通过重新设置学习速率和初始数值进行调整;或者自行判断是否线性可分;或者采用其他的分类器进行学习。
        我会在我的Github-Statistic-Learning-Method上维护扩展《统计学习方法》这本书中的其他算法,也会在我的博客ProfSnail首页继续扩展这本书中的其他内容。欢迎大家共同讨论学习。

#ifndef _PERCEPTRON
#define _PERCEPTRON
#include <cstdio>
#include <cstring>
#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;
#endif // !_PERCEPTRON

#pragma once

struct Data{
public:
	vector<double> x; // feature vector;
	int y; // instance label. 
	Data(){}
	Data(vector<double> _x, int _y): x(_x), y(_y){}
};

class Perceptron {
private:
	vector<Data> data; // train data;
	int N, n; // N: number of input data; n: vector size of each data;
	bool printLog;
	vector<double> w; // weight vector; 
	vector<double> alpha; // alpha vector: 对偶形式下的修正次数。
	vector<vector<double> > Gram;
	double b; // bias;

public:
	Perceptron(vector<Data> _data) :data(_data), N(_data.size()), n(_data[0].x.size()), w(n, 0), b(0), printLog(0), Gram(N) {}
	Perceptron(vector<Data> _data, vector<double> _w) :data(_data), N(_data.size()), n(_data[0].x.size()), w(_w), b(0), printLog(0), Gram(N) {}
	Perceptron(vector<Data> _data, vector<double> _w, double _b) :data(_data), N(_data.size()), n(_data[0].x.size()), w(_w), b(_b), printLog(0), Gram(N) {}
	Perceptron(vector<Data> _data, double _b) :data(_data), N(_data.size()), n(_data[0].x.size()), w(n, 0), b(_b), printLog(0), Gram(N) {}
	void set_Print_Log(bool _printLog) {
		printLog = _printLog;
	}
	double dotMulti(vector<double>& a, vector<double>& b) {
		// 两个向量的点乘(内积)
		double sum = 0;
		int nn = a.size();
		for (int i = 0; i < nn; i++) {
			sum += a[i] * b[i];
		}
		return sum;
	}
	double distance(int i) {
		// calculate distance from data[i] to hyperplane.
		double sum = b;
		for (int j = 0; j < n; j++) {
			sum += w[j] * data[i].x[j];
		}
		return sum;
	}

	int signx(double x) {
		// sign function.  
		if (x >= 0)return 1;
		else return -1;
	}
	
	int judgeWrong(int i, double& sum) {
		// judge whether each point got wrong classification. 
		sum = distance(i) * data[i].y;
		if (sum <= 0)return 1;
		else return 0;
	}

	double loss_hyper_distance() {
		double loss = 0;
		for (int i = 0; i < N; i++) {
			double eachLoss = 0;
			if (judgeWrong(i, eachLoss)) {
				loss += -eachLoss;
			}
		}
		return loss;
	}

	int justify_Ori_Form(double lt, int trainLimit) {
		// lt: learning rate; 0 < lt <= 1;
		int flag = 1, wrongNum=0;
		double eachLoss = 0;
		while (flag) {
			flag = 0;
			for (int i = 0; i < N; i++) {
				if (judgeWrong(i, eachLoss)) {
					for (int j = 0; j < n; j++) {
						w[j] += lt * data[i].y * data[i].x[j];
					}
					b += lt * data[i].y;
					if (printLog) {
						printf("WrongNum = %d. Data[%d] judge Wrong. w: ", ++wrongNum, i+1);
						for (int j = 0; j < n; j++) {
							printf("%.1lf ", w[j]);
						}
						printf("b: %.1lf\n", b);
						printf("w·x + b: ");
						for (int j = 0; j < n; j++) {
							printf("%.1lfx^{(%d)}+", w[j], j + 1);
						}
						printf("%.1lf\n", b);
					}
					flag = 1;
				}
			}
			/*
			一般来说,只要点集是线性可分的,就一定可以在有限次数内通过上面的随机梯度算法产生一个超平面,使得正类和
			负类完全分开。判断点集的收敛性需要用到正集和负集,两个向量组是否存线性相关,计算量较大,暂时不做预处理
			讨论,仅仅通过设置迭代次数上限进行迭代次数的限制。如果在迭代次数上限内,超平面没有完成收敛,则简单断定
			为收敛失败,检查学习速率、收敛上限或点集是否具有线性可分性。
			*/
			if (!--trainLimit) {
				return 0;
			}
		}
		if (printLog) {
			printf("Justify to accurate. w: ");
			for (int j = 0; j < n; j++) {
				printf("%.1lf ", w[j]);
			}
			printf("b: %.1lf\n", b);
		}
		printf("w·x + b: \n$");
		for (int j = 0; j < n; j++) {
			printf("%.1lfx^{(%d)}+", w[j], j + 1);
		}
		printf("%.1lf$\n", b);

		return 1;
	}

	int justify_Dual_Form(double lt, int trainLimit) {
		for(int i = 0; i < N; i++)Gram[i].resize(N);
		alpha.resize(N);
		for (int i = 0; i < N; i++) {
			for (int j = i; j < N; j++) {
				Gram[i][j] = Gram[j][i] = dotMulti(data[i].x, data[j].x);
			}
		}
		int flag = 1, wrongNum = 0;
		while (flag) {
			flag = 0;
			for (int i = 0; i < N; i++) {
				double eachCheack = b;
				for (int j = 0; j < N; j++) {
					eachCheack += alpha[j] * data[j].y * Gram[j][i];
				}
				eachCheack *= data[i].y;
				if (eachCheack <= 0) {
					alpha[i] += lt;
					b += lt * data[i].y;
					flag = 1;
					if (printLog) {
						printf("WrongNum = %d. Data[%d] judge Wrong. alpha: ", ++wrongNum, i+1);
						for (int j = 0; j < N; j++) {
							printf("%.0lf ", alpha[j]);
						}
						printf("b: %.0lf\n", b);
					}

				}
			}
			/*
			同上。
			一般来说,只要点集是线性可分的,就一定可以在有限次数内通过上面的随机梯度算法产生一个超平面,使得正类和
			负类完全分开。判断点集的收敛性需要用到正集和负集,两个向量组是否存线性相关,计算量较大,暂时不做预处理
			讨论,仅仅通过设置迭代次数上限进行迭代次数的限制。如果在迭代次数上限内,超平面没有完成收敛,则简单断定
			为收敛失败,检查学习速率、收敛上限或点集是否具有线性可分性。
			*/
			if (!--trainLimit) {
				return 0;
			}
		}
		if (printLog) {
			printf("Justify to accurate. alpha: ");
			for (int j = 0; j < N; j++) {
				printf("%.0lf ", alpha[j]);
			}
			printf("b: %.0lf\n", b);
		}
		for (int i = 0; i < N; i++) {
			for (int j = 0; j < n; j++) {
				w[j] += alpha[i] * data[i].y * data[i].x[j];
			}
		}
		printf("w·x + b: \n$");
		for (int j = 0; j < n; j++) {
			printf("%.1lfx^{(%d)}+", w[j], j + 1);
		}
		printf("%.1lf$\n", b);

		return 1;
	}

};

vector<Data> dt;
int initialData() {
	vector<double> vt(2);
	// 线性可分数据。
	
	vt[0] = 3;  vt[1] = 3;
	dt.push_back(Data(vt, 1));
	vt[0] = 4;  vt[1] = 3;
	dt.push_back(Data(vt, 1));
	vt[0] = 1;  vt[1] = 1;
	dt.push_back(Data(vt, -1));
	
	// 异或数据。线性不可分。
	/*
	vt[0] = 0;  vt[1] = 0;
	dt.push_back(Data(vt, -1));
	vt[0] = 1;  vt[1] = 1;
	dt.push_back(Data(vt, -1));
	vt[0] = 1;  vt[1] = 0;
	dt.push_back(Data(vt, 1));
	vt[0] = 0;  vt[1] = 1;
	dt.push_back(Data(vt, 1));
	*/
	return 0;
}

int main()
{
	initialData();
	Perceptron p = Perceptron(dt);
	p.set_Print_Log(1);
	// 调用原始形式的算法
	p.justify_Ori_Form(1, 100);
	// 调用对偶问题的算法。
	//p.justify_Dual_Form(1, 100);
	// 由于初始化的问题,调用一个算法之后会改变内部的超平面位置。所以只能二选一:
	// 在main函数中只能选择调用原始形式或者对偶形式两种方法之一。
	return 0;
}

        写在后面:Github上的代码已经更新,但是博文的编写工作远远比我想象中要艰难,简简单单的一篇博文,竟然写了一个星期之久,不禁对《统计学习方法》的作者李洋老师产生了由衷的敬意。如有代码需要,请移步我的github: ch02-perceptron。 欢迎大家点赞留言收藏,一键三连那!

评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

ProfSnail

谢谢老哥嗷

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值