梯度下降法及一元線性迴歸的python實現

一、梯度下降法形象解釋

　　設想我們處在一座山的半山腰的位置，現在我們需要找到一條最快的下山路徑，請問應該怎麼走？根據生活經驗，我們會用一種十分貪心的策略，即在現在所處的位置上找到一個能夠保證我們下山最快的方向，然後向著該方向行走；每到一個新位置，重複地應用上述貪心策略，我們就可以順利到達山底了。其實梯度下降法的執行過程和上述下山的例子沒有什麼區別，不同的是我們人類可以憑藉我們的感官直覺，根據所處的位置來選擇最佳的行走方向，而梯度下降法所依據的是嚴格的數學法則來進行每一步的更新。本文不再對該演算法進行嚴格的數理討論，只介紹梯度下降法進行資料擬合的流程和利用梯度下降法解決一元線性迴歸的python實現。

二、梯度下降法演算法應用流程

　　假設有一組資料X=[x₁,x₂,x₃,...]，Y=[y₁,y₂,y₃,...]，現求由X到Y的函式關係：

　　1、為所需要擬合的資料，構造合適的假設函式：y=f(x;θ)，以θ=[θ₁,θ₂,θ₃,...]為引數；

　　2、選擇合適的損失函式：cost(θ)，用損失函式來衡量假設函式對資料的擬合程度；

　　3、設定梯度下降法的學習率 α，引數的優化初始值及迭代終止條件；

　　4、迭代更新θ，直到滿足迭代終止條件，更新公式為：

　　　　θ₁=θ₁-α*dcost(θ)/dθ₁，

　　　　θ₂=θ₂-α*dcost(θ)/dθ₂，...

三、一元線性迴歸的python實現

　　下面以一個一元線性迴歸的例子來更進一步理解梯度下降法的過程。筆者通過在函式y=3*x+2的基礎之上新增一些服從均勻分佈的隨機數來構造如下的待擬合數據：X,Y，訓練資料影象如下圖1所示。假設函式為一元線性函式: y=f(x;θ,k)=θ*x+k，損失函式為:cost(θ,k)=1/2*∑(f(x_i;θ,k)-y_i),x_i屬於X，y_i屬於Y，損失函式的影象如下圖2所示。應用梯度下降法進行引數更新的過程如圖3中的藍色圓點所示。　

（1）　

（2）

（3）

　　程式原始碼如下：

 1 import numpy as np
 2 import matplotlib.pyplot as plt
 3 from mpl_toolkits.mplot3d import Axes3D
 4 
 5 np.random.seed(1)
 6 #生成樣本資料
 7 x=np.arange(-1,1,step=0.04)#自變數
 8 noise=np.random.uniform(low=-0.5,high=0.5,size=50)#噪聲
 9 y=x*3+2+noise#因變數
10 #顯示待擬合數據
11 plt.figure(1)
12 plt.xlabel('x')
13 plt.ylabel('y')
14 plt.scatter(x,y)
15 
16 #假設函式為一元線性函式：y=theta*x+k，需要求解的引數為theta和k
17 #損失函式為
18 def cost(theta, k, x, y):
19     return 1/2*np.mean((theta*x+k-y)**2)
20 
21 def cost_mesh(theta_m, k_m, x, y):
22     z_m=np.zeros((theta_m.shape[0],theta_m.shape[1]))
23     for i in range(theta_m.shape[0]):
24         for j in range(theta_m.shape[1]):
25             z_m[i,j]=cost(theta_m[i,j], k_m[i,j],x,y)
26     return z_m
27 #視覺化損失函式
28 theta_axis=np.linspace(start=0, stop=5,num=50)
29 k_axis=np.linspace(start=0, stop=5,num=50)
30 (theta_m, k_m)=np.meshgrid(theta_axis,k_axis)#網格化
31 z_m=cost_mesh(theta_m, k_m, x, y)
32 #繪製損失函式的3D影象
33 fig=plt.figure(2)
34 ax=Axes3D(fig)#為figure新增3D座標軸
35 ax.set_xlabel('theta')
36 ax.set_ylabel('k')
37 ax.set_zlabel('cost')
38 ax.plot_surface(theta_m, k_m, z_m,rstride=1, cstride=1,cmap=plt.cm.hot, alpha=0.5)#繪製3D的表面, rstide為行跨度，cstride為列跨度
39 
40 
41 #梯度下降法
42 #引數設定
43 lr=0.01#學習率
44 epoches=600#迭代次數,即迭代終止條件
45 
46 #引數初始數值
47 theta=0
48 k=0
49 
50 #迭代更新引數
51 for i in range(epoches):
52     theta_gra=np.mean((theta*x+k-y)*x)#theta梯度
53     k_gra=np.mean(theta*x+k-y)#k梯度
54     #更新梯度
55     theta-=theta_gra*lr
56     k-=k_gra*lr
57     #繪製當前引數所在的位置
58     if i%50==0:
59         ax.scatter3D(theta, k, cost(theta, k, x,y), marker='o', s=30, c='b')
60 print('最終的結果為：theta=%f, k=%f'%(theta, k))
61 plt.show()

&n