【統計學習方法】第二章 感知機

感知機模型定位：感知機屬于 二分類模型/線性模型/非概率模型/判別模型
回顧：統計學習三要素：模型+策略+演算法

演算法原理

模型

• 輸入空間/特征空間： X ? R n X \subseteq R^n X?Rn
• 輸出空間： y ∈ y \in y∈ {-1,+1}
• 輸入到輸出的映射： y = s g n ( w x + b ) y=sgn(wx+b) y=sgn(wx+b) 【sgn為符號函式】
• 假設空間：{f|f(x)=wx+b}

幾何解釋：wx+b=0是特征空間中的一個超平面S，w是該平面的法向量，b是截距；
前提假設：當資料集線性可分時，感知機才具有可用性；

策略

感知機的損失函式為誤分類的點x到超平面S的距離： 1 ∣ ∣ w ∣ ∣ ∣ w x + b ∣ \frac{1}{||w||}|wx+b| ∣∣w∣∣1?∣wx+b∣ （點到平面的距離公式），但這種含有絕對值的形式并不利于求導，因此，需要想辦法去掉絕對值；

對于誤分類的點 x i x_i xi?而言，滿足以下式子： ? y i ( w ? x i + b ) > 0 -y_i(w·x_i+b)>0 ?yi?(w?xi?+b)>0，于是，感知機的損失函式為： ? 1 ∣ ∣ w ∣ ∣ y i ( w x i + b ) -\frac{1}{||w||}y_i(wx_i+b) ?∣∣w∣∣1?yi?(wxi?+b)；

不考慮||w||，于是，就得到了感知機的風險/目標函式： L ( w , b ) = ? ∑ i y i ( w x i + b ) L(w,b)=-\sum_i y_i(wx_i+b) L(w,b)=?∑i?yi?(wxi?+b)，注意，這里的風險函式并沒有像均方誤差那樣取平均【模型的目標函式是需要根據模型的特點設定的】

演算法

感知機采用隨機梯度下降演算法進行最優解的求解；

原始形式

對L(w,b)求偏導，得到梯度：
? w L ( w , b ) = ? ∑ i y i x i \nabla_wL(w,b)=-\sum_i y_ix_i ?w?L(w,b)=?∑i?yi?xi?
? b L ( w , b ) = ? ∑ i y i \nabla_bL(w,b)=-\sum_i y_i ?b?L(w,b)=?∑i?yi?

于是，隨機選取一個誤分類點xi，w和b的更新如下：【 η \eta η為學習率】
w = w + η y i x i w=w+\eta y_ix_i w=w+ηyi?xi?
b = b + η y i b=b+\eta y_i b=b+ηyi?

對偶形式【值得仔細理解】

考慮感知機的引數更新程序，假設共進行了k次更新， k = ∑ i k i k=\sum_ik_i k=∑i?ki?，其中， k i k_i ki?為第i個點的更新次數，那么最后得到的w其實等于 w = ∑ i = 1 m α i k i y i x i w=\sum_{i=1}^m\alpha_i^{k_i}y_ix_i w=∑i=1m?αiki??yi?xi?，其中， α k i \alpha^{k_i} αki? 為對第i個樣本點的 k i k_i ki?次更新之后的引數；

直觀理解就是，對每個樣本點的更新體現在 α k i \alpha^{k_i} αki?上，而所有更新之后的樣本點之和就是w，

所以，感知機模型可定義為 y = s g n ( ∑ i = 1 m α i y i x i ? x + b ) y=sgn(\sum_{i=1}^m\alpha_iy_ix_i·x+b) y=sgn(∑i=1m?αi?yi?xi??x+b)，這里 α i \alpha_i αi?表示模型訓練后得到的最優引數

因此，我們可以將對w的更新轉換為對 α \alpha α的更新，且對誤分類點xi而言，引數更新公式為 α i = α i + η \alpha_i=\alpha_i+\eta αi?=αi?+η

注意：

• 這里的 α \alpha α是m維向量，m為輸入樣本的個數，也就是，對每個樣本，都會有一個相應的引數！
• 直觀理解引數 α \alpha α的更新：若第i個樣本被誤分類 n i n_i ni?次，則 α i \alpha_i αi?就被更新 n i n_i ni?次，每次更新，都增加 η \eta η，最后，第i個樣本對引數的貢獻為 w i = α i x i y i w_i=\alpha_ix_iy_i wi?=αi?xi?yi?，將所有樣本的引數貢獻求和，就得到了最后的w；
• 對偶形式的好處：每次進行引數更新時，無需將樣本點納入計算；

演算法收斂性——Novikoff定理

暫略

Python實作

原始形式

相關說明：

• 輸入X：m*n的矩陣，m為樣本個數，n為特征個數
• 輸出y：m*1的向量
• 引數w：n*1的向量
• 偏置b：實數

特別注意：

• 矩陣運算的實作：誰乘以誰，點乘還是矩陣乘
• 雖然說每次的引數更新是隨機選取一個誤分類點進行更新，但實際實作程序中，在一輪訓練里，一次性更新所有被誤分類的點；

'''
Author : Superpig99
Date : 2021/12/05
'''
import numpy as np

class perceptron:
    def __init__(self,learning_rate,max_epoch):
        self.lr = learning_rate # 學習率
        self.me = max_epoch # 最大的訓練次數
    # 給定X，預測y
    def predict(self,X):
        y = X @ self.w + self.b # @：矩陣乘法，維數：(m*n) * (n*1) = m*1
        y = np.where(y>=0,1,-1) # 符號函式
        return y

    def fit(self,X,y): # X是m*n的矩陣，y為m*1的向量，m為樣本個數，n為特征個數
        m,n = X.shape[0],X.shape[1]
        # 初始化
        self.w = np.zeros((n,1)) # 引數w是n*1的向量
        self.b = np.zeros(1)
        for i in range(self.me): # 開始訓練
            yhat = self.predict(X)
            wrong_index = np.where((y - yhat)!=0,1,0) # 指示矩陣，指示哪些地方預測錯了
            self.w = self.w + (self.lr*(wrong_index*y).T @ X).T # 修正w，w = w + lr * y * X，這一步很重要！值得理解
            self.b = self.b + self.lr * wrong_index.T @ y # 修正b，b = b + lr * y
            # print('epoch:',i)
            # print(self.w.T,'\n',wrong_index.T)
            print('Epoch: %d, Wrong points: %d, Error Rate: %.2f'%(i,np.sum(wrong_index),np.sum(wrong_index)/m))
            if np.sum(wrong_index)==0: # 如果全部預測正確，則訓練結束
                break
        return
    
    def evaluation(self,Yhat,Ytrue):
        if Yhat.shape == Ytrue.shape:
            acu = np.sum(np.where((Yhat - Ytrue)==0,1,0))/Ytrue.shape[0]
            return acu
        else:
            print('the shape of Yhat and Ytrue is different')
            

if __name__=='__main__':
    X = np.array([[3,3],[4,3],[1,1]])
    y = np.array([[1],[1],[-1]])
    per = perceptron(learning_rate=1,max_epoch=20)
    per.fit(X,y)
    yhat = per.predict(X)
    acu = per.evaluation(yhat,y)
    print('Accuarcy is %.2f'%acu)

重點說明：

• self.w = self.w + (self.lr*(wrong_index*y).T @ X).T該步驟含義：
  • wrong_index * y：wrong_index和y的點積（元素積），得到的是m*1的向量，含義為那些被錯誤分類的點的y值向量；
  • (wrong_index*y).T @ X)：y與X的內積，得到的是1*n的向量，含義為該輪訓練中，所有被誤分類的點的內積之和；
  • (self.lr*(wrong_index*y).T @ X).T：乘以學習率后轉置，就是該輪訓練中，w需要更新的增量；
• self.b = self.b + self.lr * wrong_index.T @ y：類推w的更新，很好理解；

對偶形式

相關說明：

• 輸入X：m*n的矩陣，m為樣本個數，n為特征個數
• 輸出y：m*1的向量
• 引數a：m*1的向量，即 α \alpha α
• 偏置b：實數

'''
Author : Superpig99
Date : 2021/12/05
'''
import numpy as np

class DaulPerceptron:
    def __init__(self,learning_rate,max_epoch):
        self.lr = learning_rate # 學習率
        self.me = max_epoch # 最大的訓練次數
    # 給定X，預測y
    def predict(self,X):
        m = X.shape[0]
        y = self.Gram @ self.c + self.b # 重點！
        y = np.where(y>=0,1,-1)
        return y

    def fit(self,X,y): # X是m*n的矩陣，y為m*1的向量，m為樣本個數，n為特征個數
        m,n = X.shape[0],X.shape[1]
        self.a = np.zeros((m,1)) # 引數a是m*1的向量
        self.b = np.zeros(1)
        self.Gram = [[0]*m for _ in range(m)] # 計算好Gram矩陣，以便以后使用
        for i in range(m):
            self.Gram[i][i] = X[i] @ X[i].T
            for j in range(i+1,m):
                self.Gram[i][j] = X[i] @ X[j].T
                self.Gram[j][i] = X[i] @ X[j].T
        for i in range(self.me): # 開始訓練
            self.c = self.a * y # 這個self.c也很重要
            yhat = self.predict(X)
            wrong_index = np.where((y - yhat)!=0,1,0) # 指示矩陣，指示哪些地方預測錯了
            self.a = self.a + self.lr*wrong_index # 修正a，a = a + lr
            self.b = self.b + self.lr*np.sum(wrong_index*y) # 修正b，b = b + lr * y
            # print('epoch:',i)
            # print(self.a.T,'\n',wrong_index.T)
            print('Epoch: %d, Wrong points: %d, Error Rate: %.2f'%(i,np.sum(wrong_index),np.sum(wrong_index)/m))
            if np.sum(wrong_index)==0: # 如果全部預測正確，則訓練結束
                break
        return
    
    def evaluation(self,Yhat,Ytrue):
        if Yhat.shape == Ytrue.shape:
            acu = np.sum(np.where((Yhat - Ytrue)==0,1,0))/Ytrue.shape[0]
            return acu
        else:
            print('the shape of Yhat and Ytrue is different')
            

if __name__=='__main__':
    X = np.array([[3,3],[4,3],[1,1]])
    y = np.array([[1],[1],[-1]])
    per = DaulPerceptron(learning_rate=1,max_epoch=20)
    per.fit(X,y)
    yhat = per.predict(X)
    acu = per.evaluation(yhat,y)
    print('Accuarcy is %.2f'%acu)

重點說明：
之前提到說，對偶形式的感知機可以寫成 y = s g n ( ∑ i = 1 m α i y i x i ? x + b ) y=sgn(\sum_{i=1}^m\alpha_iy_ix_i·x+b) y=sgn(∑i=1m?αi?yi?xi??x+b)，把式子拆看來看，這個運算式其實包含了一個Gram矩陣，元素為(xi,xj)【第i個特征向量與第j個特征向量的內積】，所以在預測的時候，計算運算式其實為y = self.Gram @ self.c + self.b，其中，self.c = self.a * y，self.c需要隨著self.a的更新而更新，這一步理解好，剩下的就都不是問題了，

總結

演算法看起來很簡單，但實作起來會發現有很多知識點會理解出錯，比如：

• 對偶形式中的引數alpha，并不是想當然的n*1的向量，而是和樣本數對應的；
• Gram矩陣是怎么來的，為什么會想到用Gram矩陣來運算，也很巧妙；

在引數更新這里，雖然表達上是說，隨機選取一個樣本點進行更新，但實際操作是每輪訓練，對所有誤分類點都進行的方法【我有看到利用for回圈對所有誤分類點進行更新的做法，但矩陣運算其實會更快】

疑問：
《統計學習方法》教材說滿足 y i ( w x i + b ) ≤ 0 y_i(wx_i+b)\leq0 yi?(wxi?+b)≤0的點都是誤分類點，教材中舉的例子也是按照 y i ( w x i + b ) ≤ 0 y_i(wx_i+b)\leq0 yi?(wxi?+b)≤0這個標準來判斷誤分類點的，但我在代碼的程序是按照預測值是否等于實際值來判斷的，所以相同的資料和初始引數下，模型更新的程序存在不同，我的疑問在于，為什么滿足等于0的點也屬于誤分類點？