[監督式]GradientDescent介绍

以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）:(真实值 - 预测值)²：

我们希望损失函数越小越好，表示预测值（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资料的平均值会变为0、 标准差变为1，能加快训练速度，促进算法的收敛。

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

特征缩放前后差异：

---

---

regularization

使多次函數降低複雜度，避免overfitting。
regularization(Lasso回歸)：

regularization(Ridge回歸、脊回歸、嶺回歸)：

• 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后可以修正过大的ｌｒ(学习率)，使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與更新參數(先求出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：迭代，每更新一次，就是一次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之间取２０个值的阵列，类型转为type， type可以为int64,int32,float64…等等。

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

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

---

---

# 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()