[ML, Python] Gradient Descent Algorithm (revision 2)

因為最近有在看Coursera複習ML, 就順便結合我這邊的一些其他資料做個整理

Gradient Descent 一個在ML中相當常用來找minimum的演算法

在進入Gradient Descent之前, 要先談一下Linear regression with one variable

Linear regression, 基本上高中數學課本就已經有教過了

給你一堆點, 要你找一條線去fit他, 可以用數學表示成 y = ax + b

那我們炫一點, 用ML常見的方式表示

x就是那些點 (input), 那 h_θ(x) 就是那條可能的線


θ其實就是以前國中高中常看到的 y = ax + b 的 a 跟 b

簡單的可能如下

Data Point:

1.1 3.17
2 3.81 
2.2 4.15 
3.1 3.75 
3.5 4.15 
4.1 4.25 
4.7 4.36 
5.3 4.56 
5.8 5.06 
6.4 5.36 

那我要用h(x)去fit他, 則, 最大的問題就是θ要怎麼找?

簡單的想法是, 我找到的f(x)要跟真實的結果, 差距最小

例如有一條線 h(x) = 1 + 2x (θ₀ = 1, θ₁ = 2)

第一個點帶進去h(1.1) = 3.2, 跟實際的3.17差了0.03, 很接近嘛

第二個點帶進去h(2) = 5, 跟實際的3.81差了1.19...後面就越差越遠了...

所以這條看來不太好, 那在繼續找其他θ來用...(找到天荒地老QO)

從這動作看出來, 如果讓計算差距產生的值越小, 就表示我的方程式越好!!

所以以數學的角度上來講要怎麼表示這件事情?

我們希望下列可以能小盡量小

看起來很可怕, 公式很噁心, 其實還好

summation裡面的就是剛剛我們做的, 把點帶入, 減去真正的y

那我們要看差異不會只看一個點, 所以才會有一個summation, 從i = 1 ~ m

意思就是總共m筆資料, 每個差距我都要加起來比較

那有個比較麻煩的問題是為甚麼會有個2次方

那簡單來說, 數學上要作最佳化, 2次方比1次方好做, 而且實作上, 2次方的最佳化結果也比較容易
(當然原因還有很多, 還有點連不連續問題...etc)

那其實這不是最完整的, 真正完整Linear Regression minimize是

甚麼!! 哪來的1 / 2m ???

這故事有點麻煩, 簡單講, 要找最小值, 那個1 / 2m不會有影響

那, 沒影響幹嘛不拿掉?

其實這是數學家的陰謀, 要讓你看不懂

其實這個數學式要作最佳化的時候, 為了讓最後的式子漂亮跟好看, 所以在這邊加的1 / 2m

也就是在推導最佳化的時候, 那個1/2會剛好不見!!!!!

那我不是教授, 所以我就不再繼續多講

簡單說, 這個看起來有點冗的數學式, 就是傳說中常見的cost function

也就是你的方程式, 拿去給cost function計算, 算出來的就是你這方程式應用在某些data上會造成的cost

這邊稍微偷懶一下, 用個代號J(θ)去表示cost function, 下面就是linear regression的cost function

(注意: 並不是cost function都長這樣, 這只是linear regression的cost設計成這樣)

---

Gradient Descent

那接下來稍微跳tone一下, 假設我們今天有個cost function (不一定是linear regression)

他有些變數, 我們希望找到一些變數, 可以讓這個function值越小越好

我們想要做的事情就是, 要不停的找θ, 並且讓J(θ)能夠越小越好

流程:

1. 初始話θ (通常都先從0開始)

2. 開始找θ, 而且希望會reduce J(θ), 希望可以找到夠小的值

那你可以無腦的窮舉的找, 但是有生之年應該找不完

我們下面簡單用個圖說明可能的情況 (隨便的function)


圖拿 matplotlib 網站上的的教學圖修改剛好可以拿來解釋

簡單講, 這是由兩個θ變數組成的所有可能的 cost 平面
(深褐色是高點, 深藍色是低點, 其餘在中間)

當然這只是部份縮圖, 這個畫面把可能的cost都畫出來了(圖裡面的數字為cost)

那, 如果把所有可能都loop過一次, 或許可以找到,  但是可能會很沒效率

假設今天我預設的點在黑黑的X那邊, 我要如何往下找到深藍色, 而不是浪費時間往深褐色的地方跑

所以我們理所當然的是要往下找, 而且找越抖的越好, 數學上要怎麼做才可以達到此效果?

那在以前的微積分(危機分!!)課程裡, 有談到要往梯度最大的地方跑, 就是要給他來做個微分阿(極值, 斜率...etc), 而且因為參數不只一種, 要多方嘗試

根據上述, 那Gradient Descent Algorithm 就是解決這些問題的演算法

PS: 不是只有Gradient Descent可以解決這個問題, 還有另外一個稱為Normal Equation

他就是一種根據J(θ)的結果, 幫你尋找適當的θ的演算法, 可以讓J(θ)達到夠小 (local minimum)

只要做以下的事情

重複下列更新 θ 直到converge

看起來好像很可怕, 其實很單純, 每次找θ的時候, 會根據J(θ)的斜率調整, 每次調整的量是α

白話一點, 就是根據J(θ)找哪個移動方向可以讓J(θ)縮小(斜率) , 然後根據移動方向移動一定距離(α), 然後在調整計算, 在看看還差目標多遠, 在繼續用同一種方式不停的, 找斜率, 移動, 找斜率, 移動...loop

Coursera上面有特別提醒一件事情, 在實作的時候, 要注意, 更新這兩個要同時更新

意思就是 θ₀, θ₁要一起更新

也就是不可以用已經更新好的θ₀去計算新的θ₁, 這樣會有問題

那回來討論一下為甚麼數學上是這樣調整

假設θ現在的值太大 (黑點)

所需要做的事情是, 往Cost斜率向下的方向調整, 那基本上就是公式的偏微部份

因為這斜率是正的, 也就是θ會往更小的地方嘗試 (公式: 正 - 正)

而要往下調整多少, 就是那個α的大小決定的, 可以稱呼他為learning rate

同理, 如果今天θ值太小

則Cost斜率向下調整, 這斜率會是負的, 所以θ會往更大的地方嘗試(公式: 正 - 負)

很直觀的方式了解這公式的作用

那這邊有個小問題, α 要怎麼取? 因為 α 就表示 θ 一次要跳多遠

很尷尬的是, 不能太大也不能太小

下面就是跳太小, 會降得很緩慢, 可能會做了老半天, 結果沒調整多少

而如果跳太大, 可能會降不了, 反而從converge變成diverse, 因為往上炸開了(跳過頭了)

對了, 稍微說明一下, Gradient Descent一般的情況下, 是只能找到Local Minimum

要符合特殊情況下, 才有可能有Global Minimum

那, 有找尋θ的方法, 目前的問題是要解linear regression, 要怎麼用?

其實也已經很單純了, 把linear regression的cost function套用Gradient Descent即可

這是微分完的結果, 另外一個θ₁在偏微的時候會多微個x_i出來

所以只要靠這個這個結果更新θ即可, 阿只要找到θ, 就可以去fit你想要fit的點囉~ (誤差可能性可以達到最小)

Coursera上補充了一點, 很有趣的是, linear regression的cost function, 在設計上, 會造成ball shape function, 也稱convex function, 會讓本來local minimum的gradient descent在這個case下, 變成global minimum (原因請看optimization class的convex相關內容)

所以linear regression的找θ的其中一個方法, 就這麼完成了!!

另外一個常見的解決方法是Normal Equation, 簡單講就是有公式使用, 不過要解反矩陣, 而且可能會有不能解的問題, 後面會有範例

Normal Equation wiki: http://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)

Implementation

首先測試data, 我拿UCI data repository的abalone來用(沒有原因, 只是隨便拿)

取第六feature為第一個feature, 跟第七個feature為label value

Download partial abalone artificial data set: dataset

我以下的程式都會利用python的numpy來使用, 這樣矩陣運算比較方便

畫圖則是會用python的matplotlib來畫圖

順邊放個原本data的分布圖

首先開頭要先設定一下

from numpy import *
import matplotlib.pyplot as plt

1. Load Dataset

首先要先讀取資料, 並且放到X matrix(training set), 以及y matrix(label set), 以方便後續使用

那讀取的function很簡單, 使用genfromtxt即可, 然後設定一下分隔是逗號( , )

要取第一維為X, 取法也很容易 data[:, 0] 就是取第一個column,

所以第二維為y, 就是用data[:, 1], 用size function來取得現在要判斷的y有幾個

data = matrix(genfromtxt('aba67', dtype=float, delimiter=','))
X = data[:, 0] 
y = data[:, 1]
m = size(y)

2. Compute Cost

再來要知道怎麼計算Cost

看圖就知道, m 是數量, h(x)是設定的θ function 然後要去扣掉 y label並且平方最後加總

說很簡單, 做起來其實也很簡單, 但是這邊有一個小tricky的東西要記住

由於 h_θ(x) = θ₀+θ₁x 是長這副模樣(迷: 有很醜嗎?), 也就是以前的y = ax + b的模樣

但是我們剛剛只有讀入X, 一個維度, 那只能搭配到θ₁, 那θ₀怎麼辦?

所以為了配合他, 我們的X要多一個維度, 全部都存放 1 !!

也就是如果原本是 matrix([[2], [3]]), 要變成 matrix([[ 1.,  2.], [ 1.,  3.]])

這樣做矩陣運算的時候才會是正確的版本

然後把重新建立好的 oneX, 跟 y , 跟 θ 一起帶入computeCost去計算即可

根據上述公式, 加上因為已經是矩陣的格式

所以基本上只要一行程式碼就可以把Cost計算給實做出來了!! 如下

def computeCost(X, y, theta):
    m = size(y);
    return (1.0 / (2.0 * m)) * sum(power(X * theta - y, 2));

oneX = concatenate((ones((m,1)), data[:,0]), axis = 1)
theta = zeros((2,1))

print computeCost(oneX, y, theta)

這邊結果應該是大約: 0.0169

3. Gradient Descent

那到最後重頭戲了, 會算Cost, 那剩下就是算Gradient了

要做Gradient, 要設定好, 還有iteration的次數, 然後呼叫做GradientDescent

根據公式, 當然是要跑 for - loop

以及重點, θ要分開更新, 所以建議是先用tmp變數去儲存修改後的θ

然後再一次更新 θ 上去, 以避免重複更新的狀況, transpose就是把 matrix 給轉置

    
def gradientDescent(X, y, theta, alpha, num_iters):
    m = size(y)
    n = size(theta)
    tmp = zeros((n,1))
    for iter in xrange(1, num_iters + 1):
        for initer in xrange(0, n):
            tmp[initer] = theta[initer] - alpha * (1.0 / m) * sum(transpose(X * theta - y) * X[:,initer])
        theta = tmp
    return theta

iterations = 3000
alpha = 0.1

theta = gradientDescent(oneX, y, theta, alpha, iterations)

將結合上述過程, 就是一個完整的Linear Regression one variable in Gradient Descent 的過程

那最後我有用plot去把最後的成果給畫出來 (藍色就是我找到的線)

以這個Dataset的角度來講, 算是滿剛好的!!

完整程式碼:

from numpy import *
import matplotlib.pyplot as plt

def computeCost(X, y, theta):
    m = size(y);
    return (1.0 / (2.0 * m)) * sum(power(X * theta - y, 2));

def gradientDescent(X, y, theta, alpha, num_iters):
    m = size(y)
    n = size(theta)
    tmp = zeros((n,1))
    for iter in xrange(1, num_iters + 1):
        for initer in xrange(0, n):
            tmp[initer] = theta[initer] - alpha * (1.0 / m) * sum(transpose(X * theta - y) * X[:,initer])
        theta = tmp
    return theta        

if __name__ == '__main__':
    data = matrix(genfromtxt('aba67', dtype=float, delimiter=','))

    X = data[:, 0] 
    y = data[:, 1]
    m = size(y)

    oneX = concatenate((ones((m,1)), data[:,0]), axis = 1)
    theta = zeros((2,1))

    print computeCost(oneX, y, theta)

    iterations = 3000
    alpha = 0.1

    theta = gradientDescent(oneX, y, theta, alpha, iterations)

    plt.plot(X, y, 'rx')
    plt.plot(oneX[:,1], oneX * theta, '-')
    plt.xlabel('Shucked weight')
    plt.ylabel('Viscera weight')
    plt.show()

上述的程式其實只要經過小修改(例如 θ 的size, 可能要做一點feature scaling...etc )

馬上就可以變成解Linear Regression multi - variable in Gradient Descent

那當然還有Normal Equation的解法

其實這超簡單實作, 所以才沒有花甚麼篇幅講

只要把最後找 θ 的方法, 改成下列即可, (公式: (X^TX)^-1X^Ty)

theta = linalg.pinv(transpose(oneX)*oneX)*transpose(oneX)*y;

最後畫出來的線跟Gradient Descent差不多

當然也是因為找到的 θ 值真的差不多

θ in Gradient Descent: [-0.00212767, 0.56028958]

θ in Normal Equation: [-0.00240189, 0.56121222]

不過那是因為這個範例簡單, 加上維度只有一個feature, 所以其實用GD的意義不大

但是如果維度一大, NE的缺點就會慢慢暴露出來!!