【Python】淺談梯度下降與實作(上):為什麼是你

設計機器學習的模型時,各種梯度下降的方法總會在我們身邊圍繞。
這次就花點時間,來仔細看看到底梯度下降都在忙些什麼。

========================================

注意事項:
   1. 為了順利顯示數學式,請先安裝 google 套件:
TeX All the Things
   
2. 本文中有些許高中數學,不會太硬,請放心咀嚼。

========================================

目次

【Python】淺談梯度下降與實作(上):為什麼是你
【Python】淺談梯度下降與實作(中):初階的變形者們
【Python】淺談梯度下降與實作(下):猙獰的變形者們

一、為什麼是你

尋找最佳解的方法有很多,相信讀者也多少聽過幾種:
像是梯度下降法、牛頓法、基因演算法、粒子演算法等等。

若搜尋梯度下降的相關教學,最常見的就是以二次函數為例來說明。
因為二次函數簡單,大家都學過,用來說明最適合不過。

若您也是剛接觸這領域的新鮮人,
可能會和之前年幼的夏恩有著相同的問題:

明明求解一次微分等於零,就能立刻知道最大值或最小值的所在位置,
為什麼要捨近求遠,還使用梯度下降法一步一步的算呢?

例如:解 $f(x)=2x^{2}-3x+1$ 的最小值,
就直接令微分等於零,求解 $f'(x)=4x-3=0$,
解得 $x=\frac{3}{4}$; 有最小值。

夏恩你看,輕鬆搞定了!
會想用梯度下降來算嗎?不會嘛!

沒錯,可以直接算當然最好,但通常算不出來。
因為實際上遇到的數學式非常難算。

假設一個基本的多變量線性回歸式:

$f(X)=w_1X_1+w_2X_2+...+w_jX_j+b ...(1)$

接著計算均方差。

$MSE=\frac{1}{n}\sum_{i=1}^{n}{(f(X^i)-Y_{gt}^i)}^2 ...(2)$

將 (1) 代入 (2) 式,以最小化均方差為目標,找到該式的最佳解。

$MSE=\frac{1}{n}\sum_{i=1}^{n}{((w_1X_1^i+w_2X_2^i+...+w_jX_j^i+b)-Y_{gt}^i)}^2 ...(3)$

恩...來算算看,看要多久才算得出來?
之後進到 deep learning 的領域,面對層層網路疊加下的產出的 loss function...
是不是覺得情不自禁地淚流滿面,不知所云?

所以區域最佳化算法的存在意義就在這裡 —
全域解計算量大,費時費力,只好以迭代區域解的方式來近似全域解

若是從全域搜索的「基因演算法」來看,找最佳解的過程可粗略分為兩個步驟:
1. 隨機給定一批權重,約莫 10~100 組(更多也可以),並計算每組權重的數值。
2. 以特定的技巧,如遺傳、交配、突變,來更新權重朝最小值靠攏。

若以「梯度下降」來說,其理念就是「走一步,看一步」。
首先隨機給一組權重,計算梯度方向,然後前進一步(即更新權重)。
接著在新的地方計算梯度方向,然後再前進一步。
反覆動作,直到收斂為止。

從上面對於兩種優化法的敘述,想必讀者也看到了些問題:
梯度下降每迭代一次,要計算 1 次全部的資料;
基因演算法每迭代一次,就要計算 100 次全部的資料!(註:取決於染色體數量)

全域最佳化算法已經很慢了,再加上資料量大的情境下,慢上加慢。

真的是算到海枯石爛,直達天荒地老!
也因此在目前主流的卷積網路優化架構中,幾乎沒有全域最佳化算法的發揮空間。

換個方面看,梯度下降法是相對簡單的 —
數學理論簡單,程式實作也簡單,計算上只算到一次微分而已。

快又有效,好東西不用嗎?

二、梯度下降(gradient descent)

梯度下降法的最基本形式,又名為:批梯度下降法(batch gradient descent, BGD)。

從初代的梯度下降法被提出後,陸續衍生許多優化算法以及各種變形體,長相也愈來愈猙獰。
既然我們想了解這傢伙是怎麼長歪的,呃不...是怎麼演化的,就得從最基本的式子看起。

梯度下降法的定義是:

$w^{j+1}=w^j-\alpha\nabla E(w),j 為迭代次數,\alpha 為學習率。$

其中 $\nabla$ 就是梯度,一次微分的意思。
若是單純的 $y=ax+b$ 的情況,一次微分就是斜率。
只是在多變量的情境中,斜率這個名稱太狹隘,於是在數學上就稱為「梯度」。

首先我們給定一個二次函數,$E(w)=w^2$。
二次函數是比較常見的誤差模型,且長得較簡單,適合說明。

在接下來的範例中,夏恩要用梯度下降來解 $E(w)$ 這個二次函數。

先複習一下,梯度的算法就是對每個變量作一次偏微分。
例如,給一函數為:$f(x,y)=x^2+2xy+y^2$。
計算其梯度為:

$\mathrm{\nabla}\ f\left(x_1,x_2\right)=\left(\frac{\partial f\left(x_1,x_2\right)}{\partial x_1},\ \frac{\partial f\left(x_1,x_2\right)}{\partial x_2}\right)=\left(2x_1+1,\ 2x_2+1\right)$

當 $(x_1, x_2)=(1, 2)$ 時,$\mathrm{\nabla}\ f\left(1,\ 2\right)=\left(3,\ 5\right)$。
上述結果用白話文來說,即 $x_1$ 方向上的斜率為 $3$;$x_2$ 方向上的斜率為 $5$ 的意思。

回過頭來看剛剛給的函數:$E(w)=w^2$,
夏恩這就來「手動」跑一次梯度下降法。

假定,第一次隨機初始化點位,$w=5$。

第一次迭代:$w$ 從 5 移動至 4

先在這邊暫停一下!
這裡有很多有趣的地方可以聊聊。

在初始化點位 (5, 25) 上,經計算後的切線斜率,也就是梯度,為 +10,這會影響幾件事情:

1. 方向:將梯度的方向取負號,就是我們想要移動的方向。

2. 大小:由於學習率固定,因此梯度值愈大,每次移動的距離愈遠!

若您已經看懂梯度所內涵的意義,那我們接著下去。

第二次迭代:$w$ 從 4 移動至 3.2

第三次迭代:$w$ 從 3.2 移動至 2.56

迭代三次,意思到了就好,畫圖也是很累人的...

這個反覆迭代的過程會一直到終止條件出現為止,例如:
1. 迭代次數達到某個值。
2. 迭代後的 loss 值間的差異小於某個數。
3. 程式執行總時間限制。

經過夏恩「手動」梯度下降,讀者們應該對於相關的流程瞭然於心了!

好吧...沒有也不勉強。

最後再來看一下學習率對於梯度下降的影響。
首先是學習率等於 0.1 的情況,​畫張圖來看看:

再來是學習率等於 0.9 的情況:

根據上圖,我們可以很明顯的看見一些現象:

學習率愈大,每次移動的距離愈長;學習率愈小則反之。
但若給值太大,則梯度下降的過程中常會有不穩定的情況發生;
若給值太小,就需要更多的迭代次數才能抵達最小值。

以下為程式碼,各位可以調整學習率來觀察不同的結果:

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

# 學習率
alpha = 0.1

# 初始位置
w_init = 5

x = np.arange(-6, 6.1, 0.1)
y = x**2
plt.plot(x, y)

plt.xlabel('w')
plt.ylabel('E(w)')
plt.title('E(w)=w^2')
plt.grid(True)

w_old = 0; w_new = 0
for i in range(1, 11):
    if i == 1:
        w_old = w_init
    else:
        w_old = w_new
    
    # (2*w_old) 為二次函數的的微分式
    w_new = w_old - alpha * (2*w_old)
    
    plt.plot((w_old, w_new), (w_old**2, w_new**2), 'ro-')
    plt.title('E(w)=w^2, error=' + '{:.2f}'.format(abs(w_new**2)))

三、實戰

明白梯度下降的基本原理後,夏恩接著就拿個實例來算一算。
資料方面,我們就先隨機產生一組數據,再來做線性回歸。

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

# 給定隨機種子
np.random.seed(1)

# =============================================== #
# functions

def getdata(n):
    # n為產生資料量

    # 給定一個固定的參數,再加上隨機變動值作為雜訊
    # 其雜訊變動值介於 -10 ~ 10 之間
    x = np.arange(-5, 5.1, 10/(n-1))
    y = 3*x + 2 + (np.random.rand(len(x))-0.5)*20
    return x, y

# =============================================== #
# 產生一組資料,共51筆。
x, y = getdata(51)

# 繪圖
plt.plot(x, y, 'bo')
plt.grid('true')
plt.ylim(-30, 30)
plt.xlabel('x')
plt.ylabel('f(x)')

好的,現在請假裝不知道剛剛設定的參數是  $f(x)=3x+2$

為了要對這組數據做線性回歸,先用我們最熟悉的方式,給定一個線性回歸式:

$f(x)=ax+b ...(4)$

上式中 $x$ 就是我們所持有的資料,$f(x)$ 是模型的預測值。
至於 $a$ 和 $b$ 則是待計算的參數,就是我們做線性回歸的目標啦!

接著設計一個計算均方差的損失函數,如此一來才能得到梯度下降的方向:

$MAE=Loss=\frac{1}{n}\sum_{i=1}^{n}{(f(x^i)-Y_t^i)}^2 ...(5),其中n為總資料筆數, Y_t 是資料的實際值。$

把上式展開,將(4)式代入(5)式,整理一下:

$Loss(a, b)=\frac{1}{n}\sum_{i=1}^{n}{((ax^i+b)-Y_t^i)}^2$

如此一來,我們就得到這次要最佳化的目標函數。

這邊您可能會問為什麼是 $f(x)-y$,而不是 $y-f(x)$ 呢?

其實計算後答案都一樣,但試想:通常$f(x)$會有很多項,$y$只有一項,
若用 $y-f(x)$,計算時就得一直變號,把所有$f(x)$內的正改負,負改正...好煩啊!
所以使用 $f(x)-y$ 就可以讓您在計算的時候省去變號的麻煩。

在開始計算前,來偷偷看一下這個損失圖長什麼樣子,
新增一個 plot_error 函數,對 $Loss(a, b)$ 作圖:

import numpy as np
import matplotlib.pyplot as plt

# 給定隨機種子
np.random.seed(1)

# =============================================== #
# functions

def getdata(n):
    # n為產生資料量

    # 給定一個固定的參數,再加上隨機變動值作為雜訊
    # 其雜訊變動值介於 -10 ~ 10 之間
    x = np.arange(-5, 5.1, 10/(n-1))
    y = 3*x + 2 + (np.random.rand(len(x))-0.5)*20
    return x, y

def plot_error(x, y):
# 繪製誤差網格圖 
   
    a = np.arange(-10, 16, 1)
    b = np.arange(-10, 16, 1)
    mesh = np.meshgrid(a, b)

    sqr_err = 0
    for xs, ys in zip(x, y):
        sqr_err += ((mesh[0]*xs + mesh[1]) - ys) ** 2
    loss = sqr_err/len(x)
    
    plt.contour(mesh[0], mesh[1], loss, 20, cmap=plt.cm.jet)
    plt.xlabel('a')
    plt.ylabel('b')
    plt.axis('scaled')
    plt.title('function loss')
    plt.colorbar()

# =============================================== #
# 產生一組資料
x, y = getdata(51)
plot_error(x, y)

上圖很明顯的是二次函數的圖形,且最小值約落在$(a,b)=(3,2)$,符合預期。

好的,言歸正傳。
還記得剛才寫的梯度下降的公式嗎?

$w^{j+1}=w^j-\alpha\nabla E(w),j 為迭代次數,\alpha 為學習率。$

現在的目標改成 $Loss(a, b)$,所以要分別計算 $a$ 和 $b$ 的梯度。

$$\mathrm{\nabla} Loss\left(a,b\right)=\left(\frac{\partial Loss\left(a,b\right)}{\partial a},\ \frac{\partial Loss\left(a,b\right)}{\partial b}\right)$$

計算後的結果如下:
需特別注意的是,由於 $Loss$ 是二次函數,計算的時候請務必要記得微積分的鏈鎖規則,

$$\frac{\partial\ Loss\left(a,\ b\right)}{\partial a}=\frac{2}{n}\sum_{i=1}^{n}{\left(ax^i+b-Y_t^i\right)\times\left(x^i\right)}$$

$$\frac{\partial\ Loss\left(a,\ b\right)}{\partial b}=\frac{2}{n}\sum_{i=1}^{n}{\left(ax^i+b-Y_t^i\right)\times\left(1\right)}$$

把上述兩式改為程式碼:

def gradient(a, b, x, y):   
    grad_a = 2 * np.mean((a*x + b - y) * (x))
    grad_b = 2 * np.mean((a*x + b - y) * (1))
    return grad_a, grad_b

有了梯度函數,我們還需要一個更新函數與均方差計算:

def update(a, b, grad_a, grad_b, alpha):
    a_new = a - alpha * grad_a
    b_new = b - alpha * grad_b
    return a_new, b_new

def mse(a, b, x, y):
    sqr_err = ((a*x + b) - y) ** 2
    loss = np.mean(sqr_err)
    return loss

這個部分,夏恩主要是想讓程式比較好懂。
若您已經看懂上述三個函數,那接著把它們結合成一個類別:

class my_BGD:    
    def __init__(self, a, b, x, y, alpha):
        self.a = a
        self.b = b
        self.x = x
        self.y = y
        self.alpha = alpha
        
        self.a_old = a
        self.b_old = b
        
        self.loss = None;
    
    # Loss function
    def mse(self):
        sqr_err = ((self.a*self.x + self.b) - self.y) ** 2
        return np.mean(sqr_err)
    
    def gradient(self):
        grad_a = 2 * np.mean((self.a*self.x + self.b - self.y) * (self.x))
        grad_b = 2 * np.mean((self.a*self.x + self.b - self.y) * (1))
        return grad_a, grad_b

    def update(self):
        # 計算梯度
        grad_a, grad_b = self.gradient()
        
        # 梯度更新
        self.a_old = self.a
        self.b_old = self.b
        self.a = self.a - self.alpha * grad_a
        self.b = self.b - self.alpha * grad_b
        self.loss = self.mse();

最後是主程式,把上面的程式碼結合起來:

# 初始設定
alpha = 0.1

# 從 -9 開始,能看得更明顯
a = -9; b = -9

# 類別 my_BGD 初始化
mlclass = my_BGD(a, b, x, y, alpha)

# 畫張圖看一下
yp = a*x + b
h, = plt.plot(x, yp, 'r-')
plt.title('Initial, loss='+'{:.2f}'.format(mlclass.mse()))

# 開始迭代 10 次
for i in range(1, 11):
    mlclass.update()
    print('iter='+str(i)+', loss='+'{:.2f}'.format(mlclass.loss))

    # 再畫張圖來看看
    yp = mlclass.a*x + mlclass.b

    h.remove()
    h, = plt.plot(x, yp, 'r-')
    plt.title('iter:'+str(i)+', loss='+'{:.2f}'.format(mlclass.loss))

看就知道學習率設太大,讀者可以自行調整學習率,大概調到 0.05 的效果就不錯。

若您也想要觀察係數$a$、$b$在誤差空間訓練的狀況,我們也可以把誤差圖畫出來看:
註:完整的程式碼請參閱附錄1。

# 初始設定
alpha = 0.1

# 從 -9 開始,能看得更明顯
a = -9; b = -9

# 初始化
mlclass = my_BGD(a, b, x, y, alpha)

plt.plot(a, b, 'ro-')
plt.title('Initial, loss='+'{:.2f}'.format(mlclass.mse())+'\na='+
          '{:.2f}'.format(a)+', b='+'{:.2f}'.format(b))

# 開始迭代
for i in range(1, 11):
    mlclass.update()
    print('iter='+str(i)+', loss='+'{:.2f}'.format(mlclass.loss))
    
    plt.plot((mlclass.a_old, mlclass.a), (mlclass.b_old, mlclass.b), 'ro-')
    plt.title('iter='+str(i)+', loss='+'{:.2f}'.format(mlclass.loss)+'\na='+
              '{:.2f}'.format(mlclass.a)+', b='+'{:.2f}'.format(mlclass.b))

在誤差空間看得比資料分布空間還要更清楚:
由於每個方向都有自己的梯度,因此兩個目標參數$a$,$b$的收斂速度不一樣。
此外學習率太大導致收斂過程大幅跳動,更甚者會直接跳過極值區域,以致無法收斂。

既然一開始就說這是「基本的」梯度下降法,能夠解決的問題也相當有限。
其本身仍有著許多非常顯而易見的缺點:

1. 當資料量過大的時候難以迭代,因為記憶體會炸裂。
2. 當愈到梯度消失的高原區,模型會找不到收斂的方向而停下來。
3. 當誤差空間非常崎嶇時,難以給定學習率。

    例如下圖是 CIFAR-10 dataset 使用 VGG-56 訓練的誤差空間:
    在這種情況下,太大的學習率就容易從山谷跳上高原;太小則是被困在區域最小值。

圖片來源:VISUALIZING THE LOSS LANDSCAPE OF NEURAL NETS

知道上述的缺點後,之後幾個變體的出現,也就理所當然了。
本文的討論先到這邊,關於梯度下降的其他變體,請於繼續下一章:

【Python】淺談梯度下降與實作(中):初階的變形者們

附錄

1. 完整的程式碼,誤差空間的收斂過程。

# -*- coding: utf-8 -*-
# 說明:誤差空間的收斂過程
#
# Shayne, 2019.09.10

import numpy as np
import matplotlib.pyplot as plt

# 給定隨機種子,使每次執行結果保持一致
np.random.seed(1)
    
def getdata(n):
    # n為產生資料量
    x = np.arange(-5, 5.1, 10/(n-1))
    # 給定一個固定的參數,再加上隨機變動值作為雜訊,其變動值介於 +-10 之間
    y = 3*x + 2 + (np.random.rand(len(x))-0.5)*20
    return x, y

def plot_error(x, y):
    a = np.arange(-10, 16, 1)
    b = np.arange(-10, 16, 1)
    mesh = np.meshgrid(a, b)

    sqr_err = 0
    for xs, ys in zip(x, y):
        sqr_err += ((mesh[0]*xs + mesh[1]) - ys) ** 2
    loss = sqr_err/len(x)
    
    plt.contour(mesh[0], mesh[1], loss, 20, cmap=plt.cm.jet)
    plt.xlabel('a')
    plt.ylabel('b')
    plt.axis('scaled')
    plt.title('function loss')

class my_BGD:    
    def __init__(self, a, b, x, y, alpha):
        self.a = a
        self.b = b
        self.x = x
        self.y = y
        self.alpha = alpha
        
        self.a_old = a
        self.b_old = b
        
        self.loss = self.mse();
    
    # Loss function
    def mse(self):
        sqr_err = ((self.a*self.x + self.b) - self.y) ** 2
        return np.mean(sqr_err)
    
    def gradient(self):
        grad_a = 2 * np.mean((self.a*self.x + self.b - self.y) * (self.x))
        grad_b = 2 * np.mean((self.a*self.x + self.b - self.y) * (1))
        return grad_a, grad_b

    def update(self):
        # 計算梯度
        grad_a, grad_b = self.gradient()
        
        # 梯度更新
        self.a_old = self.a
        self.b_old = self.b
        self.a = self.a - self.alpha * grad_a
        self.b = self.b - self.alpha * grad_b
        self.loss = self.mse();

# =============================================== #
# 隨機產生一組資料
x, y = getdata(51)

# 繪製誤差空間底圖
plot_error(x, y)

# =============================================== #

# 初始設定
alpha = 0.1

# 從 -9 開始,能看得更明顯
a = -9; b = -9

# 初始化
mlclass = my_BGD(a, b, x, y, alpha)

plt.plot(a, b, 'ro-')
plt.title('Initial, loss='+'{:.2f}'.format(mlclass.loss)+'\na='+
          '{:.2f}'.format(a)+', b='+'{:.2f}'.format(b))

# 開始迭代
for i in range(1, 11):
    mlclass.update()
    print('iter='+str(i)+', loss='+'{:.2f}'.format(mlclass.loss))
    
    plt.plot((mlclass.a_old, mlclass.a), (mlclass.b_old, mlclass.b), 'ro-')
    plt.title('iter='+str(i)+', loss='+'{:.2f}'.format(mlclass.loss)+'\na='+
              '{:.2f}'.format(mlclass.a)+', b='+'{:.2f}'.format(mlclass.b))

2. 製作 GIF 的小程式,本文的動圖是在每次迭代中存檔,最後再批次轉成 GIF 檔

# -*- coding: utf-8 -*-
# 製作 GIF

import os
import imageio

images = []
img_dir = 'path/to/your/folder'

for dirPath, _, fileNames in os.walk(img_dir):
    fileNames.sort(key=lambda x:int(x[:-4]))
    for name in fileNames:
        file_path = os.path.join(dirPath, name)
        images.append(imageio.imread(file_path))
        
imageio.mimsave('path/to/your/filename.gif', images, duration=0.5)