在求解機器學習演算法的模型引數，即無約束優化問題時，梯度下降（Gradient Descent）是最常採用的方法之一，另一種常用的方法是最小二乘法。這裡就對梯度下降法做一個完整的總結。

1. 梯度

　　　　在微積分裡面，對多元函式的引數求∂偏導數，把求得的各個引數的偏導數以向量的形式寫出來，就是梯度。比如函式f(x,y), 分別對x,y求偏導數，求得的梯度向量就是(∂f/∂x, ∂f/∂y)T,簡稱grad f(x,y)或者▽f(x,y)。對於在點(x0,y0)的具體梯度向量就是(∂f/∂x0, ∂f/∂y0)T.或者▽f(x0,y0)，如果是3個引數的向量梯度，就是(∂f/∂x, ∂f/∂y，∂f/∂z)T,以此類推。

　　　　那麼這個梯度向量求出來有什麼意義呢？他的意義從幾何意義上講，就是函式變化增加最快的地方。具體來說，對於函式f(x,y),在點(x0,y0)，沿著梯度向量的方向就是(∂f/∂x0, ∂f/∂y0)T的方向是f(x,y)增加最快的地方。或者說，沿著梯度向量的方向，更加容易找到函式的最大值。反過來說，沿著梯度向量相反的方向，也就是 -(∂f/∂x0, ∂f/∂y0)T的方向，梯度減少最快，也就是更加容易找到函式的最小值。

 　　　　

2. 梯度下降與梯度上升

　　　　在機器學習演算法中，在最小化損失函式時，可以通過梯度下降法來一步步的迭代求解，得到最小化的損失函式，和模型引數值。反過來，如果我們需要求解損失函式的最大值，這時就需要用梯度上升法來迭代了。

　　　　梯度下降法和梯度上升法是可以互相轉化的。比如我們需要求解損失函式f(θ)的最小值，這時我們需要用梯度下降法來迭代求解。但是實際上，我們可以反過來求解損失函式 -f(θ)的最大值，這時梯度上升法就派上用場了。

　　　　下面來詳細總結下梯度下降法。        

3. 梯度下降法演算法詳解

3.1 梯度下降的直觀解釋

　　　　首先來看看梯度下降的一個直觀的解釋。比如我們在一座大山上的某處位置，由於我們不知道怎麼下山，於是決定走一步算一步，也就是在每走到一個位置的時候，求解當前位置的梯度，沿著梯度的負方向，也就是當前最陡峭的位置向下走一步，然後繼續求解當前位置梯度，向這一步所在位置沿著最陡峭最易下山的位置走一步。這樣一步步的走下去，一直走到覺得我們已經到了山腳。當然這樣走下去，有可能我們不能走到山腳，而是到了某一個區域性的山峰低處。

　　　　從上面的解釋可以看出，梯度下降不一定能夠找到全域性的最優解，有可能是一個區域性最優解。當然，如果損失函式是凸函式，梯度下降法得到的解就一定是全域性最優解。

3.2 梯度下降的相關概念

　　　　在詳細瞭解梯度下降的演算法之前，我們先看看相關的一些概念。

　　　　1. 步長（Learning rate）：步長決定了在梯度下降迭代的過程中，每一步沿梯度負方向前進的長度。用上面下山的例子，步長就是在當前這一步所在位置沿著最陡峭最易下山的位置走的那一步的長度。

　　　　2.特徵（feature）：指的是樣本中輸入部分，比如2個單特徵的樣本（x(0),y(0)）,（x(1),y(1)）（x(0),y(0)）,（x(1),y(1)）,則第一個樣本特徵為x(0)x(0)，第一個樣本輸出為y(0)y(0)。

　　　　3. 假設函式（hypothesis function）：在監督學習中，為了擬合輸入樣本，而使用的假設函式，記為hθ(x)hθ(x)。比如對於單個特徵的m個樣本（x(i),y(i)）(i=1,2,...m)（x(i),y(i)）(i=1,2,...m),可以採用擬合函式如下： hθ(x)=θ0+θ1xhθ(x)=θ0+θ1x。

　　　　4. 損失函式（loss function）：為了評估模型擬合的好壞，通常用損失函式來度量擬合的程度。損失函式極小化，意味著擬合程度最好，對應的模型引數即為最優引數。線上性迴歸中，損失函式通常為樣本輸出和假設函式的差取平方。比如對於m個樣本（xi,yi）(i=1,2,...m)（xi,yi）(i=1,2,...m),採用線性迴歸，損失函式為：

             J(θ0,θ1)=∑i=1m(hθ(xi)−yi)2J(θ0,θ1)=∑i=1m(hθ(xi)−yi)2

 　　　　其中xixi表示第i個樣本特徵，yiyi表示第i個樣本對應的輸出，hθ(xi)hθ(xi)為假設函式。   

3.3 梯度下降的詳細演算法

　　　　梯度下降法的演算法可以有代數法和矩陣法（也稱向量法）兩種表示，如果對矩陣分析不熟悉，則代數法更加容易理解。不過矩陣法更加的簡潔，且由於使用了矩陣，實現邏輯更加的一目瞭然。這裡先介紹代數法，後介紹矩陣法。

 

3.3.1 梯度下降法的代數方式描述

　　　　1. 先決條件： 確認優化模型的假設函式和損失函式。

　　　　比如對於線性迴歸，假設函式表示為 hθ(x1,x2,...xn)=θ0+θ1x1+...+θnxnhθ(x1,x2,...xn)=θ0+θ1x1+...+θnxn, 其中θiθi (i = 0,1,2... n)為模型引數，xixi (i = 0,1,2... n)為每個樣本的n個特徵值。這個表示可以簡化，我們增加一個特徵x0=1x0=1 ，這樣hθ(x0,x1,...xn)=∑i=0nθixihθ(x0,x1,...xn)=∑i=0nθixi。

　　　　同樣是線性迴歸，對應於上面的假設函式，損失函式為：

           J(θ0,θ1...,θn)=12m∑j=0m(hθ(x(j)0,x(j)1,...x(j)n)−yj)2J(θ0,θ1...,θn)=12m∑j=0m(hθ(x0(j),x1(j),...xn(j))−yj)2

 

　　　　2. 演算法相關引數初始化：主要是初始化θ0,θ1...,θnθ0,θ1...,θn,演算法終止距離εε以及步長αα。在沒有任何先驗知識的時候，我喜歡將所有的θθ初始化為0， 將步長初始化為1。在調優的時候再 優化。

　　　　3. 演算法過程：

　　　　　　1）確定當前位置的損失函式的梯度，對於θiθi,其梯度表示式如下：

　　　　　　　　∂∂θiJ(θ0,θ1...,θn)∂∂θiJ(θ0,θ1...,θn)

　　　　　　2）用步長乘以損失函式的梯度，得到當前位置下降的距離，即α∂∂θiJ(θ0,θ1...,θn)α∂∂θiJ(θ0,θ1...,θn)對應於前面登山例子中的某一步。

　　　　　　3）確定是否所有的θiθi,梯度下降的距離都小於εε，如果小於εε則演算法終止，當前所有的θiθi(i=0,1,...n)即為最終結果。否則進入步驟4.

　　　　　　4）更新所有的θθ，對於θiθi，其更新表示式如下。更新完畢後繼續轉入步驟1.

　　　　　　　　θi=θi−α∂∂θiJ(θ0,θ1...,θn)θi=θi−α∂∂θiJ(θ0,θ1...,θn)

　　　　下面用線性迴歸的例子來具體描述梯度下降。假設我們的樣本是(x(0)1,x(0)2,...x(0)n,y0),(x(1)1,x(1)2,...x(1)n,y1),...(x(m)1,x(m)2,...x(m)n,ym)(x1(0),x2(0),...xn(0),y0),(x1(1),x2(1),...xn(1),y1),...(x1(m),x2(m),...xn(m),ym),損失函式如前面先決條件所述：

　　　　J(θ0,θ1...,θn)=12m∑j=0m(hθ(x(j)0,x(j)1,...x(j)n)−yj)2J(θ0,θ1...,θn)=12m∑j=0m(hθ(x0(j),x1(j),...xn(j))−yj)2。

　　　　則在演算法過程步驟1中對於θiθi 的偏導數計算如下： 　　

 　　　　∂∂θiJ(θ0,θ1...,θn)=1m∑j=0m(hθ(x(j)0,x(j)1,...x(j)n)−yj)x(j)i∂∂θiJ(θ0,θ1...,θn)=1m∑j=0m(hθ(x0(j),x1(j),...xn(j))−yj)xi(j)

　　　　由於樣本中沒有x0x0上式中令所有的xj0x0j為1.

　　　　步驟4中θiθi的更新表示式如下：

           θi=θi−α1m∑j=0m(hθ(x(j)0,x(j)1,...xjn)−yj)x(j)iθi=θi−α1m∑j=0m(hθ(x0(j),x1(j),...xnj)−yj)xi(j)

　　　　從這個例子可以看出當前點的梯度方向是由所有的樣本決定的，加1m1m 是為了好理解。由於步長也為常數，他們的乘機也為常數，所以這裡α1mα1m可以用一個常數表示。

　　　　在下面第4節會詳細講到的梯度下降法的變種，他們主要的區別就是對樣本的採用方法不同。這裡我們採用的是用所有樣本。

3.3.2 梯度下降法的矩陣方式描述

　　　　這一部分主要講解梯度下降法的矩陣方式表述，相對於3.3.1的代數法，要求有一定的矩陣分析的基礎知識，尤其是矩陣求導的知識。

　　　　1. 先決條件： 和3.3.1類似， 需要確認優化模型的假設函式和損失函式。對於線性迴歸，假設函式hθ(x1,x2,...xn)=θ0+θ1x1+...+θnxnhθ(x1,x2,...xn)=θ0+θ1x1+...+θnxn的矩陣表達方式為：

　　　　　hθ(x)=Xθhθ(x)=Xθ ，其中， 假設函式hθ(X)hθ(X)為mx1的向量,θθ為(n+1)x1的向量，裡面有n個代數法的模型引數。XX為mx(n+1)維的矩陣。m代表樣本的個數，n+1代表樣本的特徵數。

             損失函式的表示式為：J(θ)=12(Xθ−Y)T(Xθ−Y)J(θ)=12(Xθ−Y)T(Xθ−Y), 其中YY是樣本的輸出向量，維度為mx1.

　　　　2. 演算法相關引數初始化: θθ向量可以初始化為預設值，或者調優後的值。演算法終止距離εε，步長αα和3.3.1比沒有變化。

　　　　3. 演算法過程：

　　　　　　1）確定當前位置的損失函式的梯度，對於θθ向量,其梯度表示式如下：

　　　　　　　　∂∂θJ(θ)∂∂θJ(θ)

　　　　　　2）用步長乘以損失函式的梯度，得到當前位置下降的距離，即α∂∂θJ(θ)α∂∂θJ(θ)對應於前面登山例子中的某一步。

　　　　　　3）確定θθ向量裡面的每個值,梯度下降的距離都小於εε，如果小於εε則演算法終止，當前θθ向量即為最終結果。否則進入步驟4.

　　　　　　4）更新θθ向量，其更新表示式如下。更新完畢後繼續轉入步驟1.

　　　　　　　　θ=θ−α∂∂θJ(θ)θ=θ−α∂∂θJ(θ)

　　　

　　　　還是用線性迴歸的例子來描述具體的演算法過程。

　　　　損失函式對於θθ向量的偏導數計算如下：

　　　　　　∂∂θJ(θ)=XT(Xθ−Y)∂∂θJ(θ)=XT(Xθ−Y)

　　　　步驟4中θθ向量的更新表示式如下：θ=θ−αXT(Xθ−Y)θ=θ−αXT(Xθ−Y)

　　　　對於3.3.1的代數法，可以看到矩陣法要簡潔很多。這裡面用到了矩陣求導鏈式法則，和兩個矩陣求導的公式。

　　　　　　公式1：∂∂X(XXT)=2X∂∂X(XXT)=2X

　　　　　　公式2：∂∂θ(Xθ)=XT∂∂θ(Xθ)=XT

　　　　如果需要熟悉矩陣求導建議參考張賢達的《矩陣分析與應用》一書。

 

3.4 梯度下降的演算法調優

　　　　在使用梯度下降時，需要進行調優。哪些地方需要調優呢？

　　　　1. 演算法的步長選擇。在前面的演算法描述中，我提到取步長為1，但是實際上取值取決於資料樣本，可以多取一些值，從大到小，分別執行演算法，看看迭代效果，如果損失函式在變小，說明取值有效，否則要增大步長。前面說了。步長太大，會導致迭代過快，甚至有可能錯過最優解。步長太小，迭代速度太慢，很長時間演算法都不能結束。所以演算法的步長需要多次執行後才能得到一個較為優的值。

　　　　2. 演算法引數的初始值選擇。 初始值不同，獲得的最小值也有可能不同，因此梯度下降求得的只是區域性最小值；當然如果損失函式是凸函式則一定是最優解。由於有區域性最優解的風險，需要多次用不同初始值執行演算法，關鍵損失函式的最小值，選擇損失函式最小化的初值。

　　　　3.歸一化。由於樣本不同特徵的取值範圍不一樣，可能導致迭代很慢，為了減少特徵取值的影響，可以對特徵資料歸一化，也就是對於每個特徵x，求出它的期望x¯¯¯x¯和標準差std(x)，然後轉化為：

　　　　　　x−x¯¯¯std(x)x−x¯std(x)

　　　　這樣特徵的新期望為0，新方差為1，迭代次數可以大大加快。

4. 梯度下降法大家族（BGD，SGD，MBGD）

4.1 批量梯度下降法（Batch Gradient Descent）

　　　　批量梯度下降法，是梯度下降法最常用的形式，具體做法也就是在更新引數時使用所有的樣本來進行更新，這個方法對應於前面3.3.1的線性迴歸的梯度下降演算法，也就是說3.3.1的梯度下降演算法就是批量梯度下降法。　　

　　　　θi=θi−α∑j=0m(hθ(x(j)0,x(j)1,...x(j)n)−yj)x(j)iθi=θi−α∑j=0m(hθ(x0(j),x1(j),...xn(j))−yj)xi(j)

　　　　由於我們有m個樣本，這裡求梯度的時候就用了所有m個樣本的梯度資料。

4.2 隨機梯度下降法（Stochastic Gradient Descent）

　　　　隨機梯度下降法，其實和批量梯度下降法原理類似，區別在與求梯度時沒有用所有的m個樣本的資料，而是僅僅選取一個樣本j來求梯度。對應的更新公式是：

　　　　θi=θi−α(hθ(x(j)0,x(j)1,...x(j)n)−yj)x(j)iθi=θi−α(hθ(x0(j),x1(j),...xn(j))−yj)xi(j)

　　　　隨機梯度下降法，和4.1的批量梯度下降法是兩個極端，一個採用所有資料來梯度下降，一個用一個樣本來梯度下降。自然各自的優缺點都非常突出。對於訓練速度來說，隨機梯度下降法由於每次僅僅採用一個樣本來迭代，訓練速度很快，而批量梯度下降法在樣本量很大的時候，訓練速度不能讓人滿意。對於準確度來說，隨機梯度下降法用於僅僅用一個樣本決定梯度方向，導致解很有可能不是最優。對於收斂速度來說，由於隨機梯度下降法一次迭代一個樣本，導致迭代方向變化很大，不能很快的收斂到區域性最優解。

　　　　那麼，有沒有一箇中庸的辦法能夠結合兩種方法的優點呢？有！這就是4.3的小批量梯度下降法。

4.3 小批量梯度下降法（Mini-batch Gradient Descent）

　　小批量梯度下降法是批量梯度下降法和隨機梯度下降法的折衷，也就是對於m個樣本，我們採用x個樣子來迭代，1<x<m。一般可以取x=10，當然根據樣本的資料，可以調整這個x的值。對應的更新公式是：

　　　　θi=θi−α∑j=tt+x−1(hθ(x(j)0,x(j)1,...x(j)n)−yj)x(j)iθi=θi−α∑j=tt+x−1(hθ(x0(j),x1(j),...xn(j))−yj)xi(j)

5. 梯度下降法和其他無約束優化演算法的比較

　　　　在機器學習中的無約束優化演算法，除了梯度下降以外，還有前面提到的最小二乘法，此外還有牛頓法和擬牛頓法。

　　　　梯度下降法和最小二乘法相比，梯度下降法需要選擇步長，而最小二乘法不需要。梯度下降法是迭代求解，最小二乘法是計算解析解。如果樣本量不算很大，且存在解析解，最小二乘法比起梯度下降法要有優勢，計算速度很快。但是如果樣本量很大，用最小二乘法由於需要求一個超級大的逆矩陣，這時就很難或者很慢才能求解解析解了，使用迭代的梯度下降法比較有優勢。

　　　　梯度下降法和牛頓法/擬牛頓法相比，兩者都是迭代求解，不過梯度下降法是梯度求解，而牛頓法/擬牛頓法是用二階的海森矩陣的逆矩陣或偽逆矩陣求解。相對而言，使用牛頓法/擬牛頓法收斂更快。但是每次迭代的時間比梯度下降法長。

例項

Dataset
本文的資料集pga.csv包含了職業高爾夫球手的發球統計資訊，

https://raw.githubusercontent.com/huangtaosdt/Gradient-Descent-Algorithm/master/data/pga.csv

包含兩個屬性：accuracy 和 distance。accuracy 精確度描述了命中球道（ fairways hit）的比例，Distances 描述的是發球的平均距離。我們的目的是用距離來預測命中率。在高爾夫中，一個人發球越遠，那麼精度會越低。對於很多機器學習演算法來說，輸入資料前會先進行一些預處理，比如規範化，因為當計算兩個樣本的距離時，當一個屬性的取值很大，那麼這個距離會偏向取值較大的那個屬性。由於此處的精度是百分比，而距離是碼數。用到的規範化方法是每個元素減去均值然後除以標準方差。

import pandas
import matplotlib.pyplot as plt
pga = pandas.read_csv("pga.csv")

# Normalize the data  DataFrame可以用{.屬性名}來呼叫一列資料的方式，返回ndarray物件
pga.distance = (pga.distance - pga.distance.mean()) / pga.distance.std()
pga.accuracy = (pga.accuracy - pga.accuracy.mean()) / pga.accuracy.std()
print(pga.head())

plt.scatter(pga.distance, pga.accuracy)
plt.xlabel('normalized distance')
plt.ylabel('normalized accuracy')
plt.show()

#結果：
   distance  accuracy
0  0.314379 -0.707727
1  1.693777 -1.586669
2 -0.059695 -0.176699
3 -0.574047  0.372640
4  1.343083 -1.934584

Linear Model
觀察資料散點圖，發現距離和精度呈現負相關。首先用基本的線性迴歸LinearRegression 來擬合資料：

from sklearn.linear_model import LinearRegression
import numpy as np
# We can add a dimension to an array by using np.newaxis
print("Shape of the series:", pga.distance.shape) #這是一個ndarray物件，但是列數未知，因此需要人為指定一個
print("Shape with newaxis:", pga.distance[:, np.newaxis].shape) 
'''
Shape of the series: (197,)
Shape with newaxis: (197, 1)
<class 'numpy.ndarray'>
<class 'pandas.core.series.Series'>
'''
# The X variable in LinearRegression.fit() must have 2 dimensions
lm = LinearRegression()
lm.fit(pga.distance[:, np.newaxis], pga.accuracy)
theta1 = lm.coef_[0]

更簡單的方法：

# 前面之所以要新增一個維度，是因為此時屬性只有一列，但是X要求的是矩陣形（DataFrame或者二維ndarray）式，因此只要pga[['distance']]將distance放在列表中，取出來的就是一個DataFrame物件而不是Series物件。
from sklearn.linear_model import LinearRegression
import numpy as np
lm = LinearRegression()
lm.fit(pga[['distance']],pga.accuracy)
theta1 = lm.coef_[0]

Cost Function, Introduction
上述線性迴歸利用了sklearn中的包去估計模型的引數，用的是最小二乘法。最小二乘法通過矩陣計算可以很有效的擬合線性模型，並且提供確切的引數值。但是當矩陣的太大，直接用矩陣運算是不現實的，這時候就需要用一些迭代的方法來求解引數的估計值。梯度下降就是常見的一種迭代演算法。

# The cost function of a single variable linear model
def cost(theta0, theta1, x, y):
    # Initialize cost
    J = 0
    # The number of observations
    m = len(x)
    # Loop through each observation
    for i in range(m):
        # Compute the hypothesis 
        h = theta1 * x[i] + theta0
        # Add to cost
        J += (h - y[i])**2
    # Average and normalize cost
    J /= (2*m)
    return J

# The cost for theta0=0 and theta1=1
print(cost(0, 1, pga.distance, pga.accuracy))

theta0 = 100
theta1s = np.linspace(-3,2,100)
costs = []
for theta1 in theta1s:
    costs.append(cost(theta0, theta1, pga.distance, pga.accuracy))

plt.plot(theta1s, costs)

上面這段程式碼所做的工作是：將截距設定為100，係數設定為-3到2之間的100個等距離的值。然後計算每個係數所對應的模型的誤差，誤差公式如下，畫出係數與誤差的曲線圖。發現在-0.7左右，模型的誤差最小。

import numpy as np
from mpl_toolkits.mplot3d import Axes3D

theta0s = np.linspace(-2,2,100)
theta1s = np.linspace(-2,2, 100)
COST = np.empty(shape=(100,100))
# T0S：為100theta0s(theta1s的長度)行theta0s，T1S：為100列（theta0s的長度）theta1s
T0S, T1S = np.meshgrid(theta0s, theta1s)
# for each parameter combination compute the cost
for i in range(100):
    for j in range(100):
        COST[i,j] = cost(T0S[0,i], T1S[j,0], pga.distance, pga.accuracy)

# make 3d plot
fig2 = plt.figure()
ax = fig2.gca(projection='3d')
ax.plot_surface(X=T0S,Y=T1S,Z=COST)
plt.show()

上面這段程式碼首先用了一個新的函式meshgrid，引數為兩個陣列，第一個長度為m，第二個長度為n。因此返回的是第一個陣列的n行復制，以及第二個陣列的m列複製。舉個例子： 
x = [1,2,3],y=[5,6]————X=[[1,2,3],[1,2,3]],Y=[[5,5,5],[6,6,6]].然後去(X[i,j],Y[i,j])就可以得到x元素和y元素的任意組合。

上述程式碼生成了一組係數，並且將誤差與係數一起畫了一個3D圖。圖中最低的地方就是最優解。

Cost Function, Slopes
最小二乘法和梯度下降法有哪些區別？

1.本質相同：兩種方法都是在給定已知資料（independent& dependent variables）的前提下對dependent variables算出出一個一般性的估值函式。然後對給定新資料的dependent variables進行估算。 
2.目標相同：都是在已知資料的框架內，使得估算值與實際值的總平方差儘量更小（事實上未必一定要使用平方）。 
3.實現方法和結果不同：最小二乘法是直接求導找出全域性最小，是非迭代法。而梯度下降法是一種迭代法，先給定一個引數向量，然後向誤差值下降最快的方向調整引數，在若干次迭代之後找到區域性最小。梯度下降法的缺點是到最小點的時候收斂速度變慢，並且對初始點的選擇極為敏感，其改進大多是在這兩方面下功夫。

當然, 其實梯度下降法還有別的其他用處, 比如其他找極值問題. 另外, 牛頓法也是一種不錯的方法, 迭代收斂速度快於梯度下降法, 只是計算代價也比較高.

梯度下降法的步驟：初始化模型引數為w0，然後求誤差關於每個引數的偏導，偏導的反方向就是下降最快的方向，因此將引數減去偏導的值，也可以加上學習率（下降的速度）： 

計算第一個引數theta0偏導：
 

# Partial derivative of cost in terms of theta0
def partial_cost_theta0(theta0, theta1, x, y):
    # Hypothesis
    h = theta0 + theta1*x
    # Difference between hypothesis and observation
    diff = (h - y)
    # Compute partial derivative
    partial = diff.sum() / (x.shape[0])
    return partial

partial0 = partial_cost_theta0(1, 1, pga.distance, pga.accuracy)

Gradient Descent Algorithm

從之前的視覺化圖中，可以從視覺上感受到，不同的斜率和截距將帶來不同的誤差，為了減少模型的誤差，我們必須找到誤差函式中引數的最優值。下面這段程式碼就是計算梯度下降的詳細程式碼，一氣呵成，結合了前面的函式，所以在完成一個大的工程時，預先寫好一些小功能，然後在最後，拼接在一起就可以完成一個很好的功能：此處用到前面計算誤差的函式cost（），以及計算偏導的函式partial_cost_theta0（），學會這種程式設計方式。
 

# x is our feature vector -- distance
# y is our target variable -- accuracy
# alpha is the learning rate
# theta0 is the intial theta0 
# theta1 is the intial theta1
def gradient_descent(x, y, alpha=0.1, theta0=0, theta1=0):
    max_epochs = 1000 # Maximum number of iterations
    counter = 0      # Intialize a counter
    c = cost(theta1, theta0, pga.distance, pga.accuracy)  ## Initial cost
    costs = [c]     # Lets store each update
    # Set a convergence threshold to find where the cost function in minimized
    # When the difference between the previous cost and current cost 
    #        is less than this value we will say the parameters converged
    convergence_thres = 0.000001  
    cprev = c + 10   
    theta0s = [theta0]
    theta1s = [theta1]

    # When the costs converge or we hit a large number of iterations will we stop updating
    while (np.abs(cprev - c) > convergence_thres) and (counter < max_epochs):
        cprev = c
        # Alpha times the partial deriviative is our updated
        update0 = alpha * partial_cost_theta0(theta0, theta1, x, y)
        update1 = alpha * partial_cost_theta1(theta0, theta1, x, y)

        # Update theta0 and theta1 at the same time
        # We want to compute the slopes at the same set of hypothesised parameters
        #             so we update after finding the partial derivatives
        theta0 -= update0
        theta1 -= update1

        # Store thetas
        theta0s.append(theta0)
        theta1s.append(theta1)

        # Compute the new cost
        c = cost(theta0, theta1, pga.distance, pga.accuracy)

        # Store updates
        costs.append(c)
        counter += 1   # Count

    return {'theta0': theta0, 'theta1': theta1, "costs": costs}

print("Theta1 =", gradient_descent(pga.distance, pga.accuracy)['theta1'])
descend = gradient_descent(pga.distance, pga.accuracy, alpha=.01)
plt.scatter(range(len(descend["costs"])), descend["costs"])
plt.show()

以上是對資料集pga.csv操作的最終結果