Python實作多元線性回歸方程梯度下降法與求函數極值

梯度下降法

梯度下降法的基本思想可以類比為一個下山的過程。

假設這樣一個場景：一個人被困在山上，需要從山上下來(找到山的最低點，也就是山谷)。但此時山上的濃霧很大，導緻可視度很低；是以，下山的路徑就無法确定，必須利用自己周圍的資訊一步一步地找到下山的路。這個時候，便可利用梯度下降算法來幫助自己下山。怎麼做呢，首先以他目前的所處的位置為基準，尋找這個位置最陡峭的地方，然後朝着下降方向走一步，然後又繼續以目前位置為基準，再找最陡峭的地方，再走直到最後到達最低處；同理上山也是如此，隻是這時候就變成梯度上升算法了

梯度下降

梯度下降的基本過程就和下山的場景很類似。

首先，我們有一個可微分的函數。這個函數就代表着一座山。我們的目标就是找到這個函數的最小值，也就是山底。根據之前的場景假設，最快的下山的方式就是找到目前位置最陡峭的方向，然後沿着此方向向下走，對應到函數中，就是找到給定點的梯度 ，然後朝着梯度相反的方向，就能讓函數值下降的最快！因為梯度的方向就是函數之變化最快的方向(在後面會詳細解釋)

是以，我們重複利用這個方法，反複求取梯度，最後就能到達局部的最小值，這就類似于我們下山的過程。而求取梯度就确定了最陡峭的方向，也就是場景中測量方向的手段。

牛頓法

并不是所有的方程都有求根公式，或者求根公式很複雜，導緻求解困難。利用牛頓法，可以疊代求解。

原理是利用泰勒公式，在x0處展開，且展開到一階，即f(x) = f(x0)+(x－x0)f’(x0)

求解方程f(x)=0，即f(x0)+(x-x0)f’(x0)=0，求解x = x1=x0－f(x0)/f’(x0)，因為這是利用泰勒公式的一階展開，f(x) = f(x0)+(x－x0)f’(x0)處并不是完全相等，而是近似相等，這裡求得的x1并不能讓f（x）=0，隻能說f(x1)的值比f(x0)更接近f（x）=0，于是乎，疊代求解的想法就很自然了，可以進而推出x(n+1)=x(n)－f(x(n))/f’(x(n))，通過疊代，這個式子必然在f（x）=0的時候收斂。整個過程如下圖：

從本質上去看，牛頓法是二階收斂，梯度下降是一階收斂，是以牛頓法就更快。如果更通俗地說的話，比如你想找一條最短的路徑走到一個盆地的最底部，梯度下降法每次隻從你目前所處位置選一個坡度最大的方向走一步，牛頓法在選擇方向時，不僅會考慮坡度是否夠大，還會考慮你走了一步之後，坡度是否會變得更大。是以，可以說牛頓法比梯度下降法看得更遠一點，能更快地走到最底部。（牛頓法目光更加長遠，是以少走彎路；相對而言，梯度下降法隻考慮了局部的最優，沒有全局思想。）

根據wiki上的解釋，從幾何上說，牛頓法就是用一個二次曲面去拟合你目前所處位置的局部曲面，而梯度下降法是用一個平面去拟合目前的局部曲面，通常情況下，二次曲面的拟合會比平面更好，是以牛頓法選擇的下降路徑會更符合真實的最優下降路徑。

用 Python程式設計，完成函數極值的求解

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
from mpl_toolkits.mplot3d import Axes3D
import warnings

           

def f2(x1,x2):
    return x1*x1+2*x2*x2-4*x1-2*x1*x2
 
X1 = np.arange(-4,4,0.2)
X2 = np.arange(-4,4,0.2)
X1, X2 = np.meshgrid(X1, X2) # 生成xv、yv，将X1、X2變成n*m的矩陣，友善後面繪圖
Y = np.array(list(map(lambda t : f2(t[0],t[1]),zip(X1.flatten(),X2.flatten()))))
Y.shape = X1.shape # 1600的Y圖還原成原來的（40,40）
 
%matplotlib inline
#作圖
fig = plt.figure(facecolor='w')
ax = Axes3D(fig)
ax.plot_surface(X1,X2,Y,rstride=1,cstride=1,cmap=plt.cm.jet)
ax.set_title(u'$ y = x1^2+2x2^2-4x1-2x1x2 $')
plt.show()
           

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
 
def Fun(x,y):#原函數
    return x-y+2*x*x+2*x*y+y*y
 
def PxFun(x,y):#偏x導
    return 1+4*x+2*y
 
def PyFun(x,y):#偏y導
    return -1+2*x+2*y
 
#初始化
fig=plt.figure()#figure對象
ax=Axes3D(fig)#Axes3D對象
X,Y=np.mgrid[-2:2:40j,-2:2:40j]#取樣并作滿射聯合
Z=Fun(X,Y)#取樣點Z坐标打表
ax.plot_surface(X,Y,Z,rstride=1,cstride=1,cmap="rainbow")
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
 
#梯度下降
step=0.0008#下降系數
x=0
y=0#初始選取一個點
tag_x=[x]
tag_y=[y]
tag_z=[Fun(x,y)]#三個坐标分别打入表中，該表用于繪制點
new_x=x
new_y=y
Over=False
while Over==False:
    new_x-=step*PxFun(x,y)
    new_y-=step*PyFun(x,y)#分别作梯度下降
    if Fun(x,y)-Fun(new_x,new_y)<7e-9:#精度
        Over=True
    x=new_x
    y=new_y#更新舊點
    tag_x.append(x)
    tag_y.append(y)
    tag_z.append(Fun(x,y))#新點三個坐标打入表中
 
#繪制點/輸出坐标
ax.plot(tag_x,tag_y,tag_z,'r.')
plt.title('(x,y)~('+str(x)+","+str(y)+')')
plt.show()

           

用 Python程式設計，完成梯度下降法

import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
data=np.genfromtxt('D:/作業/人工實驗6.csv',delimiter=',')
x_data=data[:,:-1]
y_data=data[:,2]
#定義學習率、斜率、截據
#設方程為y=theta1x1+theta2x2+theta0
lr=0.00001
theta0=0
theta1=0
theta2=0
#定義最大疊代次數，因為梯度下降法是在不斷疊代更新k與b
epochs=10000
#定義最小二乘法函數-損失函數（代價函數）
def compute_error(theta0,theta1,theta2,x_data,y_data):
    totalerror=0
    for i in range(0,len(x_data)):#定義一共有多少樣本點
        totalerror=totalerror+(y_data[i]-(theta1*x_data[i,0]+theta2*x_data[i,1]+theta0))**2
    return totalerror/float(len(x_data))/2
#梯度下降算法求解參數
def gradient_descent_runner(x_data,y_data,theta0,theta1,theta2,lr,epochs):
    m=len(x_data)
    for i in range(epochs):
        theta0_grad=0
        theta1_grad=0
        theta2_grad=0
        for j in range(0,m):
            theta0_grad-=(1/m)*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta2)+y_data[j])
            theta1_grad-=(1/m)*x_data[j,0]*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta0)+y_data[j])
            theta2_grad-=(1/m)*x_data[j,1]*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta0)+y_data[j])
        theta0=theta0-lr*theta0_grad
        theta1=theta1-lr*theta1_grad
        theta2=theta2-lr*theta2_grad
    return theta0,theta1,theta2
#進行疊代求解
theta0,theta1,theta2=gradient_descent_runner(x_data,y_data,theta0,theta1,theta2,lr,epochs)
print('結果：疊代次數：{0} 學習率：{1}之後 a0={2},a1={3},a2={4},代價函數為{5}'.format(epochs,lr,theta0,theta1,theta2,compute_error(theta0,theta1,theta2,x_data,y_data)))
print("多元線性回歸方程為:y=",theta1,"X1+",theta2,"X2+",theta0)
#畫圖
ax=plt.figure().add_subplot(111,projection='3d')
ax.scatter(x_data[:,0],x_data[:,1],y_data,c='r',marker='o')
x0=x_data[:,0]
x1=x_data[:,1]
#生成網格矩陣
x0,x1=np.meshgrid(x0,x1)
z=theta0+theta1*x0+theta2*x1
#畫3d圖
ax.plot_surface(x0,x1,z)
ax.set_xlabel('area')
ax.set_ylabel('distance')
ax.set_zlabel("Monthly turnover")
plt.show()