美文网首页
[監督式]GradientDescent介绍

[監督式]GradientDescent介绍

作者: RJ阿杰 | 来源:发表于2018-08-21 13:15 被阅读0次

以Gradient Descent做回归分析,假设一个监督式学习的预测房价例子:
以下X为特征参数(feature):坪数,Y为标签(label):价格,图表红线为预测线。

#X = x = x1
X = [338.,333.,328.,207.,226.,25.,179.,60.,208.,606.]
#Y = y
Y = [640.,633.,619.,393.,428.,27.,193.,66.,226.,1591.]

我們假設一個function set(hypothesis set),這個function set的其中一個function 可以生成y。
我们要做Regression(回归),假设有这么一个方程(hypothesis set)代表上图的预测曲线:

single feature一次:



single feature多次:



more feature一次:

more feature多次:



补充:
Logistic Regression逻辑回归会有多个标签n(y1,y2 ... yn),同时会有多组方程,所以也会有多组b(b1,b2 ... .bn),而deep learning神经网络隐藏层也会有各自的方程,之后有空再来写Logistic Regression跟deep learning。



loss function(L)损失函数


有了这方程之后,再来我们要知道这方程预测值与真实值的误差,这时我们需要损失函数(L):(真实值 - 预测值)2

我们希望损失函数越小越好,表示预测值(z)的与真实值(y)的的差异越小,那么我们需要求得能使L最小的参数(w,b):

这个例子我们有10组,[坪数(X),价格(Y)]的资料,我们把这10组资料分别代入X,Y,然后求能使L最小的参数(w,b):

GradientDescent 梯度下降法

将X,Y代入后所得函数,以及函数曲面图:


我们要求得这曲面的最低点,我们先任意定一组初始值W0,B0,他可能会在曲面的某一点上,再来我们求这一点的梯度向量(对W,B做偏微):

我们知道斜率公式为ΔZ/Δw及ΔZ/Δb,而偏微就是点的切面斜率grad_w,grad_b:

所以grad_w为负时表示w0小,需要{w0 - grad_w}将w0调大,grad决定调整方向:

而lr为学习率,用来调整w0移动的幅度,若过大会导致无法收敛,太小会导致收敛太慢。


我们由这组梯度向量去更新w0,b0。

直到grad_w(W轴的切面斜率)与grad_b(B轴的切面斜率)为0,代表到达“局部”或“全局”最低点,此时的w0,b0就是解。

程式

再来我们以程式实作来讨论一下实作的细节,以及一些新的观念:

首先我们先导入需要用到的模组以及做资料的前处理。
matplotlib:
matplotlib是Python程式语言及其数值数学扩展包 NumPy的可视化操作界面, pylab为matplotlib基于图像处理库的接口,mpl_toolkits.mplot3d.Axes3D为matplotlib 3D图像处理库的接口。

NumPy:
NumPy是Python语言的一个扩充程式库。支援高阶大量的维度阵列与矩阵运算,此外也针对阵列运算提供大量的数学函式函式库。

SymPy:
SymPy是Python的数学符号计算库,用它可以进行数学公式的符号推导。

time:
time模块提供各种操作时间的函数。

因为我们要使用Sympy帮我们做代数运算,首先先定义变量(x,y,w,b)给sympy,定义变量的方法:



然后输入特征及标签阵列,一般会直接读取csv档(import csv、pandas、tensorflow),然后做一些资料前处理,这里先不谈。

Feature scaling(特征缩放)

再来是normalization(归一化)不同于regularization,standardization为normalization的一种,会把数据变成符合标准正态分布:
Standardization x^*= \frac{(x - mean(x))}{std(x)}


经过Standardization资料的平均值会变为0、 标准差变为1,能加快训练速度,促进算法的收敛。

訓練:x\_y\_data->norm->training->function
預測:x\_data->norm->function->norm\_y\_data
因為使用norm y訓練所以預測得到的y需要進行 reverse norm(或一開始不要對y進行normalization)。
由於這裡有對y做norm,进行预测时输入特征也需做原先归一化运算,得出的结果标签也需做反运算还原原来尺度的数据。

特征缩放前后差异:











regularization

使多次函數降低複雜度,避免overfitting。
L1 regularization(Lasso回歸):
\lambda \Sigma |w_i| + Loss\,Function
L2 regularization(Ridge回歸、脊回歸、嶺回歸):
\lambda \Sigma w_i^2 + Loss\,Function

  • underfitting、overfitting
  1. high Bias即所謂的Underfitting,因為參數過少連Training set都會有頗大的預測誤差。
  2. low Bias high Variance即所謂的Overfitting,因為參數過多導致過度符合Training set的資料特性,使得其無法預測較為普遍的資料集。


再来我们定义loss function(L)以及算出∑L,将X,Y代入L作加总。

这里使用sympy.subs(arg1,arg2)将x,y值代入,arg1为前面设置的变量(x,y),arg2为要代入的值,而subs的运算速度较差。

大量运算时改用function(f)=sympy.lambdify(arg1,arg2,arg3),arg1为变量,arg2为变量要代入的函数,arg3为要转换的代数model预设为numpy,多个变量使用f =sympy.lambdify((var1,var2),arg2,arg3)。
使用时输入f(var1,var2)将变量代入,会返回代入结果。



再来将L各自对w,b做偏微分,使用sympy.diff(function,var,order)。
然后使用numpy.random.randint(a,b,c)创建一个a~b之间,形状为c的矩阵。
ex. numpy.random.randint(1,-1,(2,2)),1 ~ -1之间,形状为2*2的矩阵。



再来就是GradientDescent实作,先设定好两个列表存更新路径,设定好要迭代的次数以及学习率,使用time计算运算时间,这里还使用了adagrad优化方法,使用adagrad后可以修正过大的lr(学习率),使GradientDescent收敛,会得到slove(w,b)。
任何需要手動調整的參數我們稱之為超參數(hyperparameters)。

adagrad优化方法

优化方法還有很多種如:Adam、Adagrad、Momentum...等等,這裡就不贅述。
下面引用几张 台大-李宏毅教授 的几张投影片来说明,详细可以去看李宏毅教授的教学课程:
https://www.youtube.com/watch?v=yKKNr-QKz2Q&index=6&list=PLJV_el3uVTsPy9oCRY30oBPNLCo89yu49

下面的g代表grad(grad_w,grad_b), η为学习率(lr),t为迭代次数。
lr(学习率)被需除以一个 σ来修正。
σ为: "计算出的grad" 与 "先前计算出的grad" 各自平方后加总,再除以加总的grad数,最后再开根号。




Batch

這是Gradient Descent求loss function常用到的方法,當data很大時我們不會一次加總所有的loss求gradient,我們會將資料batch分成幾個部分,然後在進行加總loss求出gradient與更新w,b參數(先求出gradient再加總一樣),我們這次使用的是Batch Gradient Descent所以epoch=iteration,不需另外設置epoch。

  1. Batch Gradient Descent(BGD),批梯度下降,遍歷全部數據集計算一次損失函數,進行一次參數更新,這樣得到的方向能夠更加準確的指向極值的方向,但是計算開銷大,速度慢
  2. Stochastic Gradient Descent(SGD),隨機梯度下降,對每一個樣本計算一次損失函數,進行一次參數更新,優點是速度快,缺點是方向波動大,不能準確的指向極值的方向,有時甚至兩次更新相互抵消
  3. Mini-batch Gradient Decent(MBGD),小批梯度下降,前面兩種方法的折中,把樣本數據分為若干批,裡面的資料必須是隨機的,然後分批來計算損失函數和更新參數,這樣方向比較穩定,計算開銷也相對較小

Batch Size:每一批的樣本數量,也就是每進行一次Iteration所加總樣本的的loss數量。
Iteration:迭代,w、b每更新一次,就是一次Iteration。
Epoch:樣本中的所有樣本數據被計算一次就叫做一個Epoch。



# *************************************************************
# 預測函數
def predict_fun(p_x,feature_scaling):
    if feature_scaling == 0:
        return w0 * p_x + b0
    if feature_scaling == 1:
        predict_fun_num = (p_x - np.mean(x1)) / np.std(x1)
        predict_um = (w0 *predict_fun_num+ b0) * np.std(y1) + np.mean(y1)
        return predict_um
# *************************************************************
print('\n預測800坪價格:', predict_fun(800,1))
# 更新後的w,b代入L的總loss
print('\n∑L', L_sigma_(w0,b0))
# *************************************************************

# function copy
y_copy = sp.symbols('y_copy')
x_copy = sp.symbols('x_copy')
b_copy = sp.symbols('b_copy')
w_copy = sp.symbols('w_copy')
L_copy = L_sigma
L_copy = L_copy.subs(w, w_copy).subs(b, b_copy)

W_copy = np.linspace(-3000, 3000, 20).astype(dtype='int64')
B_copy = np.linspace(-3000, 3000, 20).astype(dtype='int64')
W_copy_mesh, B_copy_mesh = np.meshgrid(W_copy, B_copy)      # w,b轉換網格矩陣
# w,b代入L求L的網格矩陣
Z = np.zeros((len(W_copy), len(B_copy)))
for mesh_I in range(len(W_copy)):
    for mesh_J in range(len(B_copy)):
        Z[mesh_J][mesh_I] = 0
        for mesh_N in range(len(X)):
            Z[mesh_J][mesh_I] = Z[mesh_J][mesh_I] + (Y[mesh_N] - B_copy[mesh_J] - W_copy[mesh_I]*X[mesh_N])**2
        Z[mesh_J][mesh_I] = Z[mesh_J][mesh_I]/len(X)

再来是做数据可视化的前处理,使用numpy.linspace(arg1,arg2,arg3).astype(dtype='type'),创建一个arg1~arg2之间取20个值的阵列,类型转为type, type可以为int64,int32,float64…等等。

使用numpy.meshgrid()建立网格矩阵:



再来求出两个网格矩阵需要对应的loss function值,最后3个网格矩阵可以用来创建3D图表。



# matplotlib plotL
# 繪製曲面圖
fig = plt.figure(figsize=(8, 4))
ax = Axes3D(fig)
plt.xlabel('w')
plt.ylabel('b')
ax.plot(W_copy, np.zeros(shape=len(W_copy)), color='k')  # 繪製3D圖的w軸
ax.plot(np.zeros(shape=len(W_copy)), B_copy, color='g')  # 繪製3D圖的b軸
ax.plot_surface(W_copy_mesh, B_copy_mesh, Z)  # 繪製曲面圖
plt.show()

# 繪製散點圖
fig1 = plt.figure(figsize=(8, 4))
ax = Axes3D(fig1)
plt.xlabel('w')
plt.ylabel('b')
ax.scatter(W_copy_mesh, B_copy_mesh, Z)  # 繪製散點圖
plt.show()

# 繪製等高線圖
fig2 = plt.figure(figsize=(8, 4))
plt.xlabel('w')
plt.ylabel('b')
plt.contourf(W_copy_mesh, B_copy_mesh, Z)                                       # 繪製等高線圖
# plt.plot(grad_path_w, grad_path_b, color='k')                                   # 繪製梯度下降w,b更新路徑直線圖
plt.scatter(grad_path_w,grad_path_b,color='r',s=4)                            # 繪製梯度下降w,b更新路徑散點圖
plt.plot([w0],[b0] , '*', ms=12, markeredgewidth=2, color='orange')       # 繪製最低點標記
plt.show()

# 繪製預測曲線
# 以原輸入的10筆X資料輸出預測對應Y值
predict_y = []
for zzz in range(len(X)):
    predict_y.append(predict_fun(X[zzz],0))
    # predict_y.append(w0 * X[zzz] + b0)
predict_y = np.array(predict_y)
# 繪製
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤(修改默認字型)
plt.rcParams['axes.unicode_minus'] = False   # 用來正常顯示負號(修改設定)
fig3 = plt.figure(figsize=(8, 4))
plt.xlabel('X:坪數')
plt.ylabel('Y:價格')
plt.scatter(X, Y, color='k')
plt.plot(X, predict_y, color='r')
plt.show()
plt.rcParams['font.sans-serif'] = ['DejaVu Sans']  # 用來正常顯示中文標籤(回復默認字型)
plt.rcParams['axes.unicode_minus'] = True  # 用來正常顯示負號(回復默認設定)

# 輸出sympy運算過程
fig4 = plt.figure(figsize=(15, 10))
plt.subplot(4, 1, 1)
plt.text(0.5, 0.5, r'L(loss function),L(w,b): $ {0} $'.format(sp.latex(L)), size=12, va='center', ha='center')
plt.subplot(4, 1, 2)
plt.text(0.5, 0.5, r'∑ L(w,b): $ {0} $'.format(sp.latex(L_sigma.expand())), size=12, va='center', ha='center')
plt.subplot(4, 1, 3)
plt.text(0.5, 0.5, r'∂L/∂w: $ {0} $'.format(sp.latex(pdw)), size=12, va='center', ha='center')
plt.subplot(4, 1, 4)
plt.text(0.5, 0.5, r'∂L/∂b: $ {0} $'.format(sp.latex(pdb)), size=12, va='center', ha='center')
plt.show()

最后使用matplotlib输出曲面图、等高线图、预测图、以及最后输出一下运算过程可视化公式。

matplotlib无法正常显示中文需要在使用到中文的图表前修改matlinplot字型设定,输出图表后在修改回来,修改后不可与特殊符号混用(中文不支持一些特殊符号)。

plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签(修改默认字型)
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号(修改设定)

输出算式使用sympy.latex可以将sympy的函数类型转成latex形式,而matlinplot可以使用自代的$$方法将latex的str字串转绘制成公式图,而这里使用python自带的str.format( 'str')将输入字串依arg顺序传入{0},{1}..。

完整code:

# 對∑L函數做偏導微(x,y未先加總)
# %matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import sympy as sp
import time
from pylab import mpl

# var
y = sp.symbols('y')
x = sp.symbols('x')
b = sp.symbols('b')
w = sp.symbols('w')

# real data X坪數 Y價格
# X = [1,9,75,83,110]
# Y = [2,18,150,170,230]
x1 = [338., 333., 328., 207., 226., 25., 179., 60., 208., 606.]
X = x1
y1 = [640., 633., 619., 393., 428., 27., 193., 66., 226., 1591.]
Y = y1

# standardization
X = (x1 - np.mean(x1)) / np.std(x1)
print('Standardization X:\n',X)
Y = (y1 - np.mean(y1)) / np.std(y1)
print('Standardization Y:\n',Y)

# lose function
L = ((y - (b + w * x)) ** 2)
print('\nL(loss function),L(w,b):\n', L)

# lose function sigma      #
sum_ = 0
for lenX0 in range(len(X)):
    reg1 = L.subs(x, X[lenX0])
    reg2 = reg1.subs(y, Y[lenX0])
    sum_ = sum_ + reg2
# ∑L:sigma L
L_sigma = sum_
L_sigma_ = sp.lambdify(( w, b), L_sigma)
print('\n∑ L(w,b):')
print('L_sigma:\n', L_sigma)
print('L_sigma_arrange:\n', L_sigma.expand())
# ---------------------   #
# 直接對L函數做偏導微 (x,y未先加總)
# ∂L/∂(w,b)
print('\nL函數做偏導微')
# ∂L/∂w
pdw = sp.diff(L, w, 1)
print('∂L/∂w:', pdw)
pdw_ = sp.lambdify((x, y, w, b), pdw)
# ∂L/∂b
pdb = sp.diff(L, b, 1)
print('∂L/∂b:', pdb)
pdb_ = sp.lambdify((x, y, w, b), pdb)

# random w0,b0
w0 = np.random.randint(-1000, 1000, 1)
w0 = w0[0]
b0 = np.random.randint(-1000, 1000, 1)
b0 = b0[0]
print('\n隨機選擇一個w,b')
print('w0:', w0)
print('b0:', b0)

# cal argmin for L
grad_path_w = [w0]
grad_path_b = [b0]
iteration = 10000
lr = 100
gradient_current_w = 0.0
gradient_current_b = 0.0
start = time.time()
for t_i in range(iteration):
    grad_w = 0
    grad_b = 0
    for batch in range(len(X)):
        grad_w = grad_w + pdw_(X[batch], Y[batch], w0, b0)                      # 由sympy計算自動代入
        grad_b = grad_b + pdb_(X[batch], Y[batch], w0, b0)                      # 由sympy計算自動代入
        # grad_w = grad_w + (-2 * X[batch] * (-b0 - w0 * X[batch] + Y[batch]))   # 由sympy計算後手動輸入
        # grad_b = grad_b + (2 * b0 + 2 * w0 * X[batch] - 2 * Y[batch])          # 由sympy計算後手動輸入
    gradient_current_w = gradient_current_w + grad_w ** 2
    gradient_current_b = gradient_current_b + grad_b ** 2
    w0 = w0 - (lr / np.sqrt(gradient_current_w)) * grad_w
    b0 = b0 - (lr / np.sqrt(gradient_current_b)) * grad_b
    grad_path_w.append(w0)
    grad_path_b.append(b0)
    # 可提早結束
    # if abs(grad_w) < 0.000005 and abs(grad_b) < 0.000005:
    #     print('break')
    #     break

print('grad_w:', grad_w)
print('grad_b:', grad_b)
print('min w:\n', w0)
print('min b:\n', b0)
end = time.time()
print('time:\n', end - start)



# *************************************************************
# 預測函數
def predict_fun(p_x,feature_scaling):
    if feature_scaling == 0:
        return w0 * p_x + b0
    if feature_scaling == 1:
        predict_fun_num = (p_x - np.mean(x1)) / np.std(x1)
        predict_um = (w0 *predict_fun_num+ b0) * np.std(y1) + np.mean(y1)
        return predict_um
# *************************************************************
print('\n預測800坪價格:', predict_fun(800,1))
# 更新後的w,b代入L的總loss
print('\n∑L', L_sigma_(w0,b0))
# *************************************************************

# function copy
y_copy = sp.symbols('y_copy')
x_copy = sp.symbols('x_copy')
b_copy = sp.symbols('b_copy')
w_copy = sp.symbols('w_copy')
L_copy = L_sigma
L_copy = L_copy.subs(w, w_copy).subs(b, b_copy)

W_copy = np.linspace(-3000, 3000, 20).astype(dtype='int64')
B_copy = np.linspace(-3000, 3000, 20).astype(dtype='int64')
W_copy_mesh, B_copy_mesh = np.meshgrid(W_copy, B_copy)      # w,b轉換網格矩陣
# w,b代入L求L的網格矩陣
Z = np.zeros((len(W_copy), len(B_copy)))
for mesh_I in range(len(W_copy)):
    for mesh_J in range(len(B_copy)):
        Z[mesh_J][mesh_I] = 0
        for mesh_N in range(len(X)):
            Z[mesh_J][mesh_I] = Z[mesh_J][mesh_I] + (Y[mesh_N] - B_copy[mesh_J] - W_copy[mesh_I]*X[mesh_N])**2
        Z[mesh_J][mesh_I] = Z[mesh_J][mesh_I]/len(X)

# matplotlib plotL
# 繪製曲面圖
fig = plt.figure(figsize=(8, 4))
ax = Axes3D(fig)
plt.xlabel('w')
plt.ylabel('b')
ax.plot(W_copy, np.zeros(shape=len(W_copy)), color='k')  # 繪製3D圖的w軸
ax.plot(np.zeros(shape=len(W_copy)), B_copy, color='g')  # 繪製3D圖的b軸
ax.plot_surface(W_copy_mesh, B_copy_mesh, Z)  # 繪製曲面圖
plt.show()

# 繪製散點圖
fig1 = plt.figure(figsize=(8, 4))
ax = Axes3D(fig1)
plt.xlabel('w')
plt.ylabel('b')
ax.scatter(W_copy_mesh, B_copy_mesh, Z)  # 繪製散點圖
plt.show()

# 繪製等高線圖
fig2 = plt.figure(figsize=(8, 4))
plt.xlabel('w')
plt.ylabel('b')
plt.contourf(W_copy_mesh, B_copy_mesh, Z)                                       # 繪製等高線圖
# plt.plot(grad_path_w, grad_path_b, color='k')                                   # 繪製梯度下降w,b更新路徑直線圖
plt.scatter(grad_path_w,grad_path_b,color='r',s=4)                            # 繪製梯度下降w,b更新路徑散點圖
plt.plot([w0],[b0] , '*', ms=12, markeredgewidth=2, color='orange')       # 繪製最低點標記
plt.show()

# 繪製預測曲線
# 以原輸入的10筆X資料輸出預測對應Y值
predict_y = []
for zzz in range(len(X)):
    predict_y.append(predict_fun(X[zzz],0))
    # predict_y.append(w0 * X[zzz] + b0)
predict_y = np.array(predict_y)
# 繪製
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用來正常顯示中文標籤(修改默認字型)
plt.rcParams['axes.unicode_minus'] = False   # 用來正常顯示負號(修改設定)
fig3 = plt.figure(figsize=(8, 4))
plt.xlabel('X:坪數')
plt.ylabel('Y:價格')
plt.scatter(X, Y, color='k')
plt.plot(X, predict_y, color='r')
plt.show()
plt.rcParams['font.sans-serif'] = ['DejaVu Sans']  # 用來正常顯示中文標籤(回復默認字型)
plt.rcParams['axes.unicode_minus'] = True  # 用來正常顯示負號(回復默認設定)

# 輸出sympy運算過程
fig4 = plt.figure(figsize=(15, 10))
plt.subplot(4, 1, 1)
plt.text(0.5, 0.5, r'L(loss function),L(w,b): $ {0} $'.format(sp.latex(L)), size=12, va='center', ha='center')
plt.subplot(4, 1, 2)
plt.text(0.5, 0.5, r'∑ L(w,b): $ {0} $'.format(sp.latex(L_sigma.expand())), size=12, va='center', ha='center')
plt.subplot(4, 1, 3)
plt.text(0.5, 0.5, r'∂L/∂w: $ {0} $'.format(sp.latex(pdw)), size=12, va='center', ha='center')
plt.subplot(4, 1, 4)
plt.text(0.5, 0.5, r'∂L/∂b: $ {0} $'.format(sp.latex(pdb)), size=12, va='center', ha='center')
plt.show()

相关文章

  • [監督式]GradientDescent介绍

    以Gradient Descent做回归分析,假设一个监督式学习的预测房价例子:以下X为特征参数(feature)...

  • [監督式、非監督]KNN、K-mean

    KNN(K Nearest Neighbor) k鄰近算法可以算是一種監督式學習算法,從部分已知的資料來推測未知的...

  • [監督式]分類

    從上篇GradientDescent延伸,這次要處理的是分類問題,從數據特徵中推斷是屬於哪一類別的標籤,例如:從身...

  • [監督式]Ensemble learning

    Ensemble learning(集成學習) 做法假設我們有多個模型(假設3個),每個模型準確率都大於50%(假...

  • 朋友,說好的日更呢

    那天朋友說,決定日更請監督。而我回覆說,寫就對了,無需監督。因為我也知道,鞭長莫及。監督有用嗎?並沒有。 這是第二...

  • [監督式]PLA(Perceptron Learning Alg

    架構圖 為群體參數,為一個函數輸入得到 (資料)為由群體取出的樣本參數加上雜訊,為一個函數,輸入得到(預測) 為推...

  • [監督式]SVM(Support Vector Machines

    SVM(Support Vector Machines) 建議可以先讀PLA(感知器) 這邊與PLA不同的是,我們...

  • [監督式]seq2seq

    預備知識RNN、LSTM之前花了一些時間讀了seq2seq+attention,發現沒做筆記過一陣子真的很容易忘,...

  • [監督式]Deep learning(深度學習)

    Deep learning Neural Network 每一層都可以使用不同的,sigmoid、ReLU、lea...

  • [監督式]NLP基礎-RNN、LSTM

    RNN RNN常用來處理序列數據,例如我們有10年每月的天氣紀錄要來預測作物收成,我們必須要知道是否有持續好幾天的...

网友评论

      本文标题:[監督式]GradientDescent介绍

      本文链接:https://www.haomeiwen.com/subject/uoudiftx.html