Linear Regression 簡介

假設我們有 \(m\) 筆資料，每筆資料都包含 \(n\) 個特徵，記為 \(x\)，每筆資料都有 1 個正確答案，記為 \(y\)。若給我們一筆新的資料，我們希望預測這筆資料的正確答案。

\(x\) 與 \(y\) 之間可能存在著某種關聯性 \(f\)，使得 \(y=f(x)\)，但我們不知道這個 \(f\) 是什麼，因此我們希望利用這些現存的資料找出一個 \(h(x)\)，使得 \(h\) 越接近 \(f\) 越好。

Model Representation

Linear Regression (線性迴歸) 就是找出一條直線，使得這條線能夠盡可能的擬合所有的訓練資料。

\[h_\theta(x)=\theta_0+\theta_1 x_1+\theta_2 x_2+...+\theta_n x_n\]

\(\theta\) 是這條直線的係數，也是我們要想辦法計算出來的，這個模型的參數。

那麼假設我們找到了一組參數，我們要如何保證在給定這組參數的情況下，模型能夠準確地用來表示這些資料呢？

Cost Function

我們希望資料透過這個模型的計算得出來的結果，與實際的值誤差越小越好，因此我們可以採用 Mean Square Error (MSE) 來計算誤差。

而這個誤差我們稱為 cost function (代價函數)，有些人稱作 loss function (損失函數) 或 objective function (目標函數)，其實指的都是同個意思。

下方的式子就是我們使用的 cost function：

\[J(\theta)=\frac{1}{2m}\sum_{i=1}^m(h_\theta(x^{(i)})-y^{(i)})^2\]

這裡有一點需要注意的是，我們多除了一個 \(2\)，這是因為待會的步驟要計算這個 cost function 的微分，微分時平方項會變成 \(\times2\)，剛好可以與多除的 \(2\) 抵銷掉，因此能夠減少計算量。

Gradient Descent

既然我們有了一個 cost function 能夠衡量一組參數的好壞，那我們剩下要做的事情就是想辦法找到一組最好的參數使得 cost function 最小。

我們當然可以用暴力法一一嘗試各種參數的組合，但是這個有點方法不切實際，因為我們不知道參數的範圍在哪裡，同時這個方法的效率也相當差。

我們可以用 Gradient Descent 來逐步更新參數：

\[\theta_j:=\theta_j-\alpha\frac{\partial}{\partial \theta_j}J(\theta)\]

\(\dfrac{\partial}{\partial \theta_j}J(\theta)\) 是 cost function 對參數 \(\theta_j\) 的偏微分，代表這一點的梯度方向，而由於梯度方向是指向高點，而我們要尋找的是 cost function 的最小值，所以我們會在這一項前面加上一個負號。

我們將這個梯度乘上一個 learning rate \(\alpha\)，這個值通常小於 1，因為我們希望每次參數更新的幅度不要太大，才能夠緩慢的收斂至極值，不過若 learning rate 過小也容易導致收斂在 local optimal。

而 gradient 的計算如下：

\[\dfrac{\partial}{\partial \theta_j}J(\theta)=\dfrac{1}{m}\sum_{i=1}^m(h_\theta(x^{(i)})-y^{(i)})x_j^{(i)}\]

Feature Scaling

若不同 feature 的範圍差距過大，在 gradient descent 的時候收斂的速度會較慢，這時我們可以將 feature 進行正規化，例如將特徵規範在 -1 到 1 之間。

而正規化的方法有很多種，這裡我們使用的是 Mean Normalization：

\[x_i=\dfrac{x_i-\mu_i}{s_i}\]

其中 \(s_i=x_{max}-x_{min}\)，結果會使得 \(-1\leq x_i\leq1\)，並且 mean = 0。

Normal Equation

除了用 gradient descent 之外，也可以使用 Normal Equation，這個方法能夠直接得出最佳的 \(\theta\) 值，這邊我們僅列出式子，不探討其背後的原理。

\[\theta=(X^TX)^{-1}X^Ty\]

實作

載入我們唯一需要用到的函式庫：numpy。

1

import numpy as np

定義一個類別及所需的變數：

1
2
3
4
5
6
7
8
9

class LinearRegression():
    def __init__(self, num_iteration=100, learning_rate=1e-1, feature_scaling=True):
        self.num_iteration = num_iteration
        self.learning_rate = learning_rate
        self.feature_scaling = feature_scaling
        self.M = 0 # normalize mean
        self.S = 1 # normalize range
        self.W = None
        self.cost_history = np.empty(num_iteration)

Gradient Descent 的方法：

為了加快運算的速度以及方便表示資料，我們會把資料向量化，也就是全部的資料都以numpy的矩陣表示。

因為 feature 的數量為 \(n\) 個，因此參數的數量會是 \(n + 1\) 個：

\[\theta=\begin{bmatrix} \theta_0 \\ \theta_1 \\ \vdots \\ \theta_n \end{bmatrix}_{(n+1)\times 1} \]

但輸入資料為一個 \(m\times n\) 的矩陣：

\[X=\begin{bmatrix} x_1^{(1)} & x_2^{(1)} & \dots & x_n^{(1)} \\ x_1^{(2)} & x_2^{(2)} & \dots & x_n^{(2)} \\ \vdots & \vdots & \vdots & \vdots \\ x_1^{(m)} & x_2^{(m)} & \dots & x_n^{(m)} \end{bmatrix}_{m\times n} \]

因此需要新增一行全部都為 \(1\) 的 feature 與 \(\theta_0\) 做計算：

\[X=\begin{bmatrix} 1 & x_1^{(1)} & x_2^{(1)} & \dots & x_n^{(1)} \\ 1 & x_1^{(2)} & x_2^{(2)} & \dots & x_n^{(2)} \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 1 & x_1^{(m)} & x_2^{(m)} & \dots & x_n^{(m)} \end{bmatrix}_{m\times (n+1)} \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

def fit(self, X, y):
    # m 為資料筆數，n 為特徵數量
    if X.ndim == 1:
        X = X.reshape(X.shape[0], 1)
    m, n = X.shape

    # 是否進行正規化
    if self.feature_scaling:
        X = self.normalize(X)

    # 在 X 左方加入一行 1 對應到參數 theta 0
    X = np.hstack((np.ones((m, 1)), X))

    y = y.reshape(y.shape[0], 1)

    self.W = np.zeros((n+1,1))

    # 每個 iteration 逐步更新參數
    for i in range(self.num_iteration):
        y_hat = X.dot(self.W)
        cost = self.cost_function(y_hat, y, m)
        self.cost_history[i] = cost
        self.gradient_descent(X, y_hat, y, m)

實作fit函式用到的方法：

1
2
3
4
5
6
7
8
9
10
11
12
13

def normalize(self, X):
    self.M = np.mean(X, axis=0)
    self.S = np.max(X, axis=0) - np.min(X, axis=0)
    return (X - self.M) / self.S

def cost_function(self, y_hat, y, m):
    return 1/(2*m) * np.sum((y_hat - y)**2)

def compute_gradient(self, X, y_hat, y, m):
    return 1/m * np.sum((y_hat - y) * X, axis=0).reshape(-1,1)

def gradient_descent(self, X, y_hat, y, m):
    self.W -= self.learning_rate * self.compute_gradient(X, y_hat, y, m)

預測：

1
2
3
4
5
6
7
8
9
10
11
12

def predict(self, X):
    if X.ndim == 1:
        X = X.reshape(X.shape[0], 1)
    m, n = X.shape

    if self.normalize:
        X = (X - self.M) / self.S

    X = np.hstack((np.ones((m, 1)), X))

    y_hat = X.dot(self.W)
    return y_hat

Normal Equation 的方法：

1
2
3
4
5
6
7
8
9
10

def normal_equation(self, X, y):
    if X.ndim == 1:
        X = X.reshape(X.shape[0], 1)
    m, n = X.shape

    X = np.hstack((np.ones((m, 1)), X))

    y = y.reshape(y.shape[0], 1)

    self.W = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)

實驗

隨機產生 30 筆單一特徵的資料：

不使用 feature scaling

使用 feature scaling

其實可以發現在 univariate 的 feature 上，不論使用 feature scaling 與否，都不會對結果影響太大。

使用 Normal Equation：

隨機產生 100 筆三個特徵的資料，前三個 column 為 feature，最後一個 column 為 label，y 的值透過三個特徵經過某些運算後加上一些隨機值得到：

1
2
3
4
5
6

54.31 3.92 -2.06 130.03
34.79 3.63 -2.49 91.56
87.91 2.99 -2.97 197.75
15.88 3.4 -2.22 52.14
96.59 2.83 -1.81 209.64
...

不進行 feature scaling，則 learning rate 必須設定小一點，否則容易造成梯度爆炸。即便沒有梯度爆炸，也通常需要更多的 iteration 才會收斂。

下圖使用的 learning rate = 1e-5。

若有進行 feature scaling，則 learning rate 可以設定大一點，收斂速度會相對較快。

下圖使用的 learning rate = 0.1。