1. 程式人生 > >Logistic迴歸之梯度上升優化演算法(三)

Logistic迴歸之梯度上升優化演算法(三)

Logistic迴歸之梯度上升優化演算法(三)

1、改進的隨機梯度上升演算法

前面兩節講了Logistic迴歸以及裡面常用的梯度上升優化演算法來找到最佳迴歸係數。但是梯度上升優化演算法的計算量很大,每次更新迴歸係數時都需要遍歷整個資料集。下面給出之前所講的梯度上升演算法:

def gradAscent(dataMatIn, classLables):
    dataMatrix = np.mat(dataMatIn)  #轉換成numpy的mat
    # print(dataMatrix)
    labelMat =  np.mat(classLables).transpose() #轉換成numpy的mat,並進行轉置
    # print(labelMat)
    m, n =np.shape(dataMatrix)#返回dataMatrix的大小。m為行,n為列
    alpha = 0.001  #移動補償,也就是學習速率,控制更新的幅度
    maxCycles = 500 #最大迭代次數
    weights = np.ones((n,1))
    # print(weights)
    for k in range(maxCycles):
        h = sigmoid(dataMatrix *weights) #梯度上升向量公式
        # print(h)
        error = labelMat - h
        weights = weights + alpha * dataMatrix.transpose()*error
    return weights.getA()  #將矩陣轉換為陣列,返回權重陣列

該演算法在對一個100個樣本的資料集處理時,dataMatrix是一個100*3的矩陣,每次計算h的時候,都要計算dataMatrix*weights這個矩陣乘法運算,要進行100*3次乘法運算和100*2次加法運算。在更新weights時也需要對整個資料集進行處理,對於大樣本是不可取的,為此使用一種隨機梯度上升演算法來簡化演算法。

1.1、隨機梯度上升演算法

我們直接給出程式碼:

def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = np.shape(dataMatrix) #返回dataMatrix的大小 m為行數 n為列數
    weights = np.ones(n) #[1. 1. 1.] 初始化引數 型別是array所以注意dataMatrix型別也應該是np.array

    for j in range(numIter):
        dataIndex = list(range(m))  # [0,1,...99]
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01 #降低alpha的大小,每次減小1/(j+i)
            randIndex = int(random.uniform(0,len(dataIndex)))#隨機選取樣本
            h = sigmoid(sum(dataMatrix[randIndex]*weights))#選擇隨機選取的一個樣本,計算h。對應位置相乘求和
            error = classLabels[randIndex] - h
            weights = weights+ alpha* error*dataMatrix[randIndex]#更新迴歸係數
            del(dataIndex[randIndex])#刪除已經使用的樣本
    return weights

該演算法主要是在以下兩方面做了改進。第一,alpha在每次迭代的時候都會調整,隨著迭代次數不斷減小,但永遠不會減小到0。這樣做的原因是為了保證在多次迭代之後新資料仍然具有一定的影響。如果要處理的問題是動態變化的,那麼可以適當加大上述常數項,來確保新的值獲得更大的迴歸係數。而且在降低alpha的函式中,alpha每次減少1/(j+i),其中j是迭代次數,i是樣本的下標。這樣當j<<max(i)時,alpha就不是嚴格下降的。避免參數的嚴格下降也常見於模擬退火演算法等其他優化演算法中。第二,這裡通過每次隨機選取一個樣本來更新迴歸係數。這種方法可以有效的減少週期性的波動。

程式碼如下:

import matplotlib.pyplot as plt
import numpy as np
import random
from matplotlib.font_manager import FontProperties
'''
函式說明:載入資料
Parameters:
    None
Returns:
    dataMat - 資料列表
    labelMat - 標籤列表
'''
def loadDataSet():
    dataMat = []  # 建立資料列表
    labelMat = []  # 建立標籤列表
    fr = open('testSet.txt')  # 開啟檔案
    for line in fr.readlines():  # 逐行讀取
        lineArr = line.strip().split()  # 去回車,放入列表
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])  # 新增資料
        labelMat.append(int(lineArr[2]))  # 新增標籤
    fr.close()  # 關閉檔案
    return dataMat, labelMat  # 返回
'''
函式說明:繪製資料集
Parameters:
    None
Returns:
    None
'''
def plotDataSet(weights):
    dataMat, labelMat = loadDataSet()  # 載入資料集
    dataArr = np.array(dataMat)  # 轉換成numpy的array陣列
    n = np.shape(dataMat)[0]  # 資料個數,即行數
    xcord1 = [] ; ycord1 = []  # 正樣本
    xcord2 = [] ; ycord2 = []  # 負樣本
    for i in range(n):
        if int(labelMat[i]) == 1: #1為正樣本
            xcord1.append(dataMat[i][1])
            ycord1.append(dataMat[i][2])
            # xcord1.append(dataArr[i, 1]);ycord1.append(dataArr[i, 2])
        else:                     #0為負樣本
            xcord2.append(dataMat[i][1])
            ycord2.append(dataMat[i][2])
            # xcord2.append(dataArr[i, 1]);ycord2.append(dataArr[i, 2])
    fig = plt.figure()
    ax = fig.add_subplot(111)   #新增subplot
    ax.scatter(xcord1,ycord1,s=20,c='red',marker = 's', alpha=.5,label ='1') #繪製正樣本
    ax.scatter(xcord2,ycord2,s=20,c='green',marker = 's', alpha=.5,label ='0') #繪製正樣本
    x = np.arange(-3.0,3.0,0.1)
    y = (-weights[0] - weights[1] * x) / weights[2]
    ax.plot(x,y)
    plt.title('DataSet') #繪製title
    plt.xlabel('x'); plt.ylabel('y') #繪製label
    plt.legend()
    plt.show()
'''
函式說明:sigmodi函式
Paremeters:
    inX - 資料
Returns:
    sigmoid函式
'''
def sigmoid(inX):
    return 1.0/(1 + np.exp(-inX))
'''
函式說明:梯度上升演算法
Parameters:
    dataMatIn - 資料集
    classLables - 資料標籤
Returns:
    weights.getA() - 求得的權重陣列(最優引數)
'''
# def gradAscent(dataMatIn, classLables):
#     dataMatrix = np.mat(dataMatIn)  #轉換成numpy的mat
#     # print(dataMatrix)
#     labelMat =  np.mat(classLables).transpose() #轉換成numpy的mat,並進行轉置
#     # print(labelMat)
#     m, n =np.shape(dataMatrix)#返回dataMatrix的大小。m為行,n為列
#     alpha = 0.001  #移動補償,也就是學習速率,控制更新的幅度
#     maxCycles = 500 #最大迭代次數
#     weights = np.ones((n,1))
#     # print(weights)
#     for k in range(maxCycles):
#         h = sigmoid(dataMatrix *weights) #梯度上升向量公式
#         # print(h)
#         error = labelMat - h
#         weights = weights + alpha * dataMatrix.transpose()*error
#     return weights.getA()  #將矩陣轉換為陣列,返回權重陣列
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = np.shape(dataMatrix) #返回dataMatrix的大小 m為行數 n為列數
    weights = np.ones(n) #[1. 1. 1.] 初始化引數 型別是array所以注意dataMatrix型別也應該是np.array
    for j in range(numIter):
        dataIndex = list(range(m))  # [0,1,...99]
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01 #降低alpha的大小,每次減小1/(j+i)
            randIndex = int(random.uniform(0,len(dataIndex)))#隨機選取樣本
            h = sigmoid(sum(dataMatrix[randIndex]*weights))#選擇隨機選取的一個樣本,計算h。對應位置相乘求和
            error = classLabels[randIndex] - h
            weights = weights+ alpha* error*dataMatrix[randIndex]#更新迴歸係數
            del(dataIndex[randIndex])#刪除已經使用的樣本
    return weights
if __name__ == '__main__':
    np.set_printoptions(suppress=True)#關閉科學技術法
    dataMat,labelMat = loadDataSet()
    weights = stocGradAscent1(np.array(dataMat),labelMat)
    plotDataSet(weights)

程式碼執行結果:

1.2、迴歸係數與迭代次數的關係

可以看出隨機梯度演算法得到的迴歸係數與初始結果所差無幾。接下來我們研究迭代次數和迴歸係數的關係。

from matplotlib.font_manager import FontProperties
import matplotlib.pyplot as plt
import numpy as np
import random
'''
函式說明:載入資料
Parameters:
    無
Returns:
    dataMat - 資料列表
    labelMat - 標籤列表
'''
def loadDataSet():
    dataMat = [] #建立資料列表
    labelMat = [] #建立標籤列表
    fr = open('testSet.txt') #開啟檔案
    for line in fr.readlines(): #逐行讀取
        lineArr = line.strip().split() #去回車,放入列表
        dataMat.append([1.0 ,float(lineArr[0]), float(lineArr[1])])
        labelMat.append(int(lineArr[2]))#新增標籤
    fr.close()
    return dataMat,labelMat
'''
函式說明:sigmoid函式
Parameters:
    inX - 資料
Returns:
    sigmoid函式
'''
def sigmoid(inX):
    return 1.0/(1+np.exp(-inX))
'''
函式說明:梯度上升演算法
Parameters:
    dataMatIn - 資料集
    classLables - 資料標籤
Returns:
    weight.getA() -  求得的權重陣列(最優引數)
    weight_array - 每次更新的迴歸係數
'''
def gradAscent(dataMatIn ,classLabels):
    dataMatrix = np.mat(dataMatIn) #轉換成numpy的mat
    labelMat = np.mat(classLabels).transpose()#轉換成numpy的mat並進行轉置
    m,n = np.shape(dataMatrix)#獲得dataMatrix的行列
    alpha = 0.01
    maxCycles = 500
    weights = np.ones((n,1))
    weights_array = np.array([])
    for k in range(maxCycles):
        h = sigmoid(dataMatrix * weights)
        error = labelMat - h
        weights = weights + alpha * dataMatrix.transpose() * error
        weights_array = np.append(weights_array,weights)#把weights插入weights_array最後面
    weights_array = weights_array.reshape(maxCycles,n)
    return weights.getA(), weights_array#將矩陣轉換為陣列並返回
'''
函式說明:改進的隨機梯度上升演算法
Parameters:
    dataMatrix - 資料陣列
    classLabels - 資料標籤
Returns:
    weights - 求得的迴歸係數陣列(最優引數)
    weights_array - 每次更新的迴歸係數
'''
def stocGradAscent1(dataMatrix ,classLabels,numIter=150):
    m,n = np.shape(dataMatrix)#返回dataMatrix的大小
    weights = np.ones(n) #引數初始化
    weights_array = np.array([]) #儲存每次更新的迴歸係數
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01#降低alpha的大小
            randIndex = int(random.uniform(0,len(dataIndex)))#隨機選取樣本
            h = sigmoid(sum(dataMatrix[randIndex]*weights))
            error = classLabels[randIndex] - h
            weights = weights + alpha * error * dataMatrix[randIndex]
            weights_array = np.append(weights_array,weights,axis = 0)
            del(dataIndex[randIndex])
    weights_array = weights_array.reshape(numIter*m,n)
    return weights,weights_array
'''
函式說明:繪製迴歸係數與迭代次數的關係
Parameters:
    weights_array1 - 迴歸係數陣列1
    weights_array2 - 迴歸係數陣列2
Returns:
    None
'''
def plotWeihts(weights_array1,weights_array2):
    #設定漢子格式
    font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=14)
    #將fig畫布分割成1行1列,不共享x軸和y軸,fig畫布的大小為(13,8)
    #講nrow=3,nclos=2時,代表fig畫布被分成6個區域,axs[0][0]表示第一行第一列
    fig,axs = plt.subplots(nrows=3,ncols=2,sharex=False,sharey=False,figsize=(20,10))
    x1 = np.arange(0,len(weights_array1))#等差數列
    #繪製w0與迭代次數的關係
    axs[0][0].plot(x1,weights_array1[:,0]) #每一行的第一個元素
    axs0_title_text = axs[0][0].set_title('梯度上升演算法:迴歸係數與迭代次數關係',FontProperties=font)
    axs0_ylabel_text = axs[0][0].set_ylabel('W0',FontProperties=font)
    plt.setp(axs0_title_text,size=20,weight='bold',color = 'black')
    plt.setp(axs0_ylabel_text,size=20,weight='bold',color = 'black')
  #繪製w1與迭代次數的關係
    axs[1][0].plot(x1,weights_array1[:,1]) #每一行的第一個元素
    axs0_ylabel_text = axs[1][0].set_ylabel('W1',FontProperties=font)
    plt.setp(axs0_ylabel_text,size=20,weight='bold',color = 'black')
  #繪製w2與迭代次數的關係
    axs[2][0].plot(x1,weights_array1[:,2]) #每一行的第一個元素
    axs0_xlabel_text = axs[2][0].set_xlabel('迭代次數',FontProperties=font)
    axs0_ylabel_text = axs[2][0].set_ylabel('W2',FontProperties=font)
    plt.setp(axs0_xlabel_text,size=20,weight='bold',color = 'black')
    plt.setp(axs0_ylabel_text,size=20,weight='bold',color = 'black')

    x2 = np.arange(0, len(weights_array2))
    # 繪製w0與迭代次數的關係
    axs[0][1].plot(x2, weights_array2[:, 0])
    axs0_title_text = axs[0][1].set_title(u'改進的隨機梯度上升演算法:迴歸係數與迭代次數關係', FontProperties=font)
    axs0_ylabel_text = axs[0][1].set_ylabel(u'W0', FontProperties=font)
    plt.setp(axs0_title_text, size=20, weight='bold', color='black')
    plt.setp(axs0_ylabel_text, size=20, weight='bold', color='black')
    # 繪製w1與迭代次數的關係
    axs[1][1].plot(x2, weights_array2[:, 1])
    axs1_ylabel_text = axs[1][1].set_ylabel(u'W1', FontProperties=font)
    plt.setp(axs1_ylabel_text, size=20, weight='bold', color='black')
    # 繪製w2與迭代次數的關係
    axs[2][1].plot(x2, weights_array2[:, 2])
    axs2_xlabel_text = axs[2][1].set_xlabel(u'迭代次數', FontProperties=font)
    axs2_ylabel_text = axs[2][1].set_ylabel(u'W2', FontProperties=font)
    plt.setp(axs2_xlabel_text, size=20, weight='bold', color='black')
    plt.setp(axs2_ylabel_text, size=20, weight='bold', color='black')
    plt.show()
if __name__ == '__main__':
    dataMat , labelMat = loadDataSet()
    weights1,weights_array1 = gradAscent(dataMat,labelMat)
    weights2, weights_array2 = stocGradAscent1(np.array(dataMat), labelMat)
    plotWeihts(weights_array1,weights_array2)




執行結果如下:

由於改進的隨機梯度上升演算法,隨機選取樣本點,所以每次的執行結果有些不同但是大體趨勢是一樣的。從圖中我們可以看出隨機梯度上升演算法的瘦臉效果更好。我們一共有100個樣本點,改進的隨機梯度上升演算法迭代次數為150.而上圖顯示150000次迭代次數的原因是,使用一次樣本就更新一下回歸係數,要迭代150次。所以150*100=150000次。迭代150次,更新1.5萬次迴歸引數。從右圖我們可以知道基本迭代2000次即迭代20次迴歸係數已經收斂了,比初始演算法收斂要更快一些。

參考文獻: