梯度下降從放棄到入門

梯度

在向量微積分中，標量場的梯度是一個向量場。標量場中某一點上的梯度指向標量場增長最快的方向，梯度的長度是這個最大的變化率。

在三維直角坐標係中表示為：

$$
\nabla \varphi = \left( \frac{\partial \varphi}{\partial x}, \frac{\partial \varphi}{\partial y}, \frac{\partial \varphi}{\partial z}\right)
= \frac{\partial \varphi}{\partial x}\vec{i} + \frac{\partial \varphi}{\partial y}\vec{j} + \frac{\partial \varphi}{\partial z}\vec{k}
$$

梯度下降

梯度下降是一個用來求函數最小值的算法，在機器學習所有優化算法當中，梯度下降應該算是最受關注並且用的最多的一個算法。

梯度下降背後的思想是：開始時我們隨機選取一個參數的組合$(\theta_0,\theta_1,\cdots,\theta_n)$，計算代價函數$J(\theta)$，然後我們尋找下一個能讓代價函數值下降最多的參數組合。我們持續這麼做，直到找到一個局部最小值，因為我們並沒有嚐試完所有的參數組合，所以不能確定我們得到局部最小值是否便是全局最小值，不同的初始參數組合，可能會找到不同的局部最小值。如果我們知道代價函數是一個凸函數，那麼我們就可以保證我們找到的局部最小點就是全局最小點。

通過梯度的定義，我們知道，梯度的方向是增長（上升）最快的方向，那麼梯度的反方向便是讓代價函數下降程度最大的方向。

定義$\alpha$為學習率（learning rate），它決定了我們沿著能讓代價函數下降程度最大的方向更新的步長。

$$ \theta_j = \theta_j - \alpha \frac{\partial J(\theta)}{\partial \theta_j}$$

每走一步，我們都要重新計算函數在當前點的梯度，然後選擇梯度的反方向作為走下去的方向。隨著每一步迭代，梯度不斷地減小，到最後減小為零。

線性回歸

假定我們有一批數據點，需要根據這些點找出最適合的線性方程$y = ax + b$。

X = [1, 2, 3, 4, 5, 6, 7, 8, 9]
Y = [0.199, 0.389, 0.580, 0.783, 0.980, 1.177, 1.380, 1.575, 1.771]

也就是要找出最合適的參數$(a, b)$，使得直線到所有點的距離盡可能小。他可以轉化為求解代價函數1

$$ J(a, b) = \sum_{i=1}^m\left(ax^{(i)} + b - y^{(i)}\right)^2 $$

的最小值。

我們也可以令$J(a, b)$為誤差平方和（SSE）的平均數，為了後續求導計算方便，我們還可以再把他除以2，即代價函數2

$$ J(a, b) = \frac{1}{2m}\sum_{i=1}^m\left(ax^{(i)} + b - y^{(i)}\right)^2 $$

如果$(a, b)$在代價函數2上取得最小值，也必在代價函數1上取的最小值。

其中$$x^{(i)}, y^{(i)}$$為已知數。

$$
\frac{\partial}{\partial a}J(a,b)
= \frac{\partial}{\partial a}\left(\frac{1}{2m}\sum_{i=1}^m\left(ax^{(i)} + b - y^{(i)}\right)^2\right)
= \frac{1}{m}\sum_{i=1}^m\left(ax^{(i)} + b - y^{(i)}\right)x^{(i)}
$$

$$
\frac{\partial}{\partial b}J(a,b)
= \frac{\partial}{\partial b}\left(\frac{1}{2m}\sum_{i=1}^m\left(ax^{(i)} + b - y^{(i)}\right)^2\right)
= \frac{1}{m}\sum_{i=1}^m\left(ax^{(i)} + b - y^{(i)}\right)
$$

利用梯度下降法，我們很容易實現對應的迭代解法。

$$ a := a - \alpha \frac{\partial}{\partial a}J(a,b) $$

$$ b := b - \alpha \frac{\partial}{\partial b}J(a,b) $$

def gd(X, Y, alpha=0.01, epsilon=1e-8):
    m = len(X)
    a, b, sse2 = 0, 0, 0
    while True:
        grad_a, grad_b = 0, 0
        for i in range(m):
            diff = a * X[i] + b - Y[i]
            grad_a += X[i] * diff
            grad_b += diff

        grad_a = grad_a / m
        grad_b = grad_b / m

        a -= alpha * grad_a
        b -= alpha * grad_b

        sse = 0
        for j in range(m):
            sse += (a  X[j] + b - Y[j])  2 / (2  m)

        if abs(sse2 - sse) < epsilon:
            break
        else:
            sse2 = sse
    return a, b

測試代碼如下:

a, b = gd(X, Y, 0.05, 1e-6)

print 'y = {0} * x + {1}'.format(a, b)
# y = 0.193958973089 * x + 0.0161212366801
x = np.array(X)
plt.plot(x, Y, 'o', label='Original data', markersize=5)
plt.plot(x, a * x + b, 'r', label='Fitted line')
plt.show()

多變量線性回歸

上例，我們演示了一元線性回歸梯度下降的迭代解法，更一般地，我們考慮n個變量的情況，我們有m條記錄。

$$ \left[
\begin{array}{ccccc|c}
x^{(1)}_{0} & x^{(1)}_{1} & x^{(1)}_{2} & \cdots & x^{(1)}_{n} & y^{(1)} \
x^{(2)}_{0} & x^{(2)}_{1} & x^{(2)}_{2} & \cdots & x^{(2)}_{n} & y^{(2)} \
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \
x^{(m)}_{0} & x^{(m)}_{1} & x^{(m)}_{2} & \cdots & x^{(m)}_{n} & y^{(m)} \
\end{array}
\right]
$$

為了不失一般性和簡化公式，我們令 $\left.x_0^{(i)} = 1\right|_{i \in (1,2,\cdots,m)}$

我們需要找到一個假設$h$，使得應用到上述數據，其代價函數最小。

$$ h_{\theta}(x) = \theta_0x_0 + \theta_1x_1 + \theta_2x_2 + \cdots + \theta_nx_n $$

或者

h_{\theta}(x) = \sum_{j=0}^n\theta_jx_j

\begin{bmatrix}
\theta_0 & \theta_1 & \cdots & \theta_n
\end{bmatrix}
\cdot
\begin{bmatrix}
x_0 \
x_1 \
\vdots \
x_n \
\end{bmatrix}
= \theta^TX
$$

其中，$x_0 = 1$, $\theta_0$為偏置。

$$ J(\theta) = \frac{1}{2m}\sum_{i=1}^m\left(h_{\theta}\left(x^{(i)}\right) - y^{(i)}\right)^2 $$

則$J(\theta)$梯度為

$$ \nabla J(\theta) =
\begin{bmatrix}
\frac{\partial J(\theta)}{\partial\theta_0} \
\frac{\partial J(\theta)}{\partial\theta_1} \
\frac{\partial J(\theta)}{\partial\theta_2} \
\vdots \
\frac{\partial J(\theta)}{\partial\theta_n} \
\end{bmatrix}
$$

$$
\frac{\partial J(\theta)}{\partial\theta_j}
= \frac{\partial}{\partial\theta_j}\left(\frac{1}{2m}\sum_{i=1}^m\left(h_{\theta}\left(x^{(i)}\right) - y^{(i)}\right)^2\right)
= \frac{1}{m}\sum_{i=1}^m\left(h_{\theta}\left(x^{(i)}\right) - y^{(i)}\right)x_j^{(i)}
$$

利用梯度下降法，我們很容易實現對應的迭代解法。

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

def gd(X, Y, alpha=0.01, epsilon=1e-6):
    m, n = len(X), len(X[0])
    theta = [1 for i in range(n)]
    sse2 = 0
    while True:
        gradient = [0 for i in range(n)]
        for j in range(n):
            for i in range(m):
                hypothesis = sum(X[i][jj] * theta[jj] for jj in range(n))
                loss = hypothesis - Y[i]
                gradient[j] += X[i][j] * loss
            gradient[j] = gradient[j] / m

        for j in range(n):
            theta[j] = theta[j] - alpha * gradient[j]

        sse = 0
        for i in range(m):
            loss = sum(X[i][jj] * theta[jj] for jj in range(n)) - Y[i]
            sse += loss ** 2
        sse = sse / (2 * m)

        if abs(sse2 - sse) < epsilon:
            break
        else:
            sse2 = sse
    return theta

之所以不在同一個循環裏麵同時計算梯度和更新$\theta$，是因為正確的做法需要保證$\theta$同時更新。

# 取x_0 = 1
X = [(1, 1.), (1, 2.), (1, 3.), (1, 4.), (1, 5.), (1, 6.), (1, 7.), (1, 8.), (1, 9.)]
Y = [0.199, 0.389, 0.580, 0.783, 0.980, 1.177, 1.380, 1.575, 1.771]

b, a = gd(X, Y, 0.05, 1e-6)

print 'y = {0} * x + {1}'.format(a, b)
# y = 0.193953964855 * x + 0.01615274985

更優雅的實現

之前一直有個疑問，不清楚為什麼用顯卡可以提高優化算法的執行效率，直到我接觸到了如下代碼：

loss = np.dot(X, theta) - Y
gradient = np.dot(X.T, loss) / m
theta -= alpha * gradient

利用矩陣運算，代碼可以變的相當簡潔，而且可以利用顯卡多核的優勢。

推導過程如下：

LOSS

\begin{bmatrix}
loss^{(1)} \
loss^{(2)} \
\vdots \
loss^{(m)} \

\end{bmatrix}

\begin{bmatrix}
h_{\theta}(x^{(1)}) - y^{(1)} \
h_{\theta}(x^{(2)}) - y^{(2)} \
\vdots \
h_{\theta}(x^{(m)}) - y^{(m)} \

\end{bmatrix}

\begin{bmatrix}
\sum_{j=0}^nx^{(1)}_{j}\theta_j \
\sum_{j=0}^nx^{(2)}_{j}\theta_j \
\vdots \
\sum_{j=0}^nx^{(m)}_{j}\theta_j \

\end{bmatrix}

\begin{bmatrix}
y^{(1)} \
y^{(2)} \
\vdots \
y^{(m)} \
\end{bmatrix}
$$

$$ = \begin{bmatrix}
x^{(1)}_{0} & x^{(1)}_{1} & x^{(1)}_{2} & \cdots & x^{(1)}_{n} \
x^{(2)}_{0} & x^{(2)}_{1} & x^{(2)}_{2} & \cdots & x^{(2)}_{n} \
\vdots & \vdots & \vdots & \ddots & \vdots \
x^{(m)}_{0} & x^{(m)}_{1} & x^{(m)}_{2} & \cdots & x^{(m)}_{n} \
\end{bmatrix}
\cdot
\begin{bmatrix}
\theta_0 \
\theta_1 \
\vdots \
\theta_n \

\end{bmatrix}

\begin{bmatrix}
y^{(1)} \
y^{(2)} \
\vdots \
y^{(m)} \
\end{bmatrix}
= X \cdot \theta - Y
$$

$$
GRADIENT =
\begin{bmatrix}
\frac{\partial}{\partial\theta_0}J(\theta) \
\frac{\partial}{\partial\theta_1}J(\theta) \
\vdots \
\frac{\partial}{\partial\theta_n}J(\theta) \

\end{bmatrix}

\begin{bmatrix}
\sum_{i=1}^m(h_\theta(x^{(i)} - y^{(i)})x^{(i)}_{0} \
\sum_{i=1}^m(h_\theta(x^{(i)} - y^{(i)})x^{(i)}_{1} \
\vdots \
\sum_{i=1}^m(h_\theta(x^{(i)} - y^{(i)})x^{(i)}_{n} \

\end{bmatrix}

\begin{bmatrix}
\sum_{i=1}^mloss^{(i)}x^{(i)}_{0} \
\sum_{i=1}^mloss^{(i)}x^{(i)}_{1} \
\vdots \
\sum_{i=1}^mloss^{(i)}x^{(i)}_{n} \
\end{bmatrix}
$$

$$
= \begin{bmatrix}
x^{(1)}_{0} & x^{(2)}_{0} & \cdots & x^{(m)}_{0} \
x^{(1)}_{1} & x^{(2)}_{1} & \cdots & x^{(m)}_{1} \
\vdots \
x^{(1)}_{n} & x^{(2)}_{n} & \cdots & x^{(m)}_{n} \
\end{bmatrix}
\cdot
\begin{bmatrix}
loss^{(1)} \
loss^{(2)} \
\vdots \
loss^{(m)} \
\end{bmatrix}
= X^T \cdot LOSS
$$

完整代碼如下

import numpy as np

def gd(X, Y, alpha=0.01, epsilon=1e-6):
    m, n = np.shape(X)
    theta = np.ones(n)
    sse2 = 0
    Xt = np.transpose(X)
    while True:
        hypothesis = np.dot(X, theta)
        loss = hypothesis - Y

        sse = np.dot(loss.T, loss) / (2 * m)
        if abs(sse2 - sse) < epsilon:
            break
        else:
            sse2 = sse

        gradient = np.dot(Xt, loss) / m
        theta -= alpha * gradient

    return theta

# 取x_0 = 1
X = [(1, 1.), (1, 2.), (1, 3.), (1, 4.), (1, 5.), (1, 6.), (1, 7.), (1, 8.), (1, 9.)]
Y = [0.199, 0.389, 0.580, 0.783, 0.980, 1.177, 1.380, 1.575, 1.771]

b, a = gd(X, Y, 0.05, 1e-6)

print 'y = {0} * x + {1}'.format(a, b)
# y = 0.193953964855 * x + 0.01615274985

最終實現

在上麵的使用上，我們需要調用方主動為每一個$x_0$賦值為1，其實這一部分完全可以在算法裏麵實現，簡化調用方的使用成本。

import numpy as np

def bgd(X, Y, alpha=0.01, epsilon=1e-8, trace=True):
    m = len(X)
    _X = np.column_stack((np.ones(m), X))
    m, n = np.shape(_X)
    theta, sse2, cnt = np.ones(n), 0, 0
    Xt = _X.T

    while True:
        loss = np.dot(_X, theta) - Y

        sse = np.dot(loss.T, loss) / (2 * m)
        if abs(sse - sse2) < epsilon:
            break
        else:
            sse2 = sse

        if trace:
            print "[ Epoch {0} ] theta = {1}, loss = {2}, error = {3})".format(cnt, theta, loss, sse)

        gradient = np.dot(Xt, loss) / m
        theta -= alpha * gradient
        cnt += 1

    return theta

X = [1, 2, 3, 4, 5, 6, 7, 8, 9]
Y = [1.99, 3.89, 5.80, 7.83, 9.80, 11.77, 13.80, 15.75, 17.71]

b, a = gd(X, Y, 0.05, 1e-6)

print 'y = {0} * x + {1}'.format(a, b)

最後更新：2017-11-06 14:05:14

梯度下降從放棄到入門

梯度

梯度下降

線性回歸

多變量線性回歸

h_{\theta}(x) = \sum_{j=0}^n\theta_jx_j

更優雅的實現

LOSS

\end{bmatrix}

\end{bmatrix}

\end{bmatrix}

\end{bmatrix}

\end{bmatrix}

\end{bmatrix}

最終實現

上一篇：量邦科技馮永昌：AI讓量化投資的戰爭升級，交易麵或許不改變其零和博弈的性質｜人工智能研習社

下一篇：深到骨子裏的自律，是每周堅持刷幾篇最新論文 | PaperDaily #10

相關內容

熱門內容

最新內容

梯度下降從放棄到入門

梯度

梯度下降

線性回歸

多變量線性回歸

h_{\theta}(x) = \sum_{j=0}^n\theta_jx_j

更優雅的實現

LOSS

\end{bmatrix}

\end{bmatrix}

\end{bmatrix}

\end{bmatrix}

\end{bmatrix}

\end{bmatrix}

最終實現

上一篇： 量邦科技馮永昌：AI讓量化投資的戰爭升級，交易麵或許不改變其零和博弈的性質｜人工智能研習社

下一篇： 深到骨子裏的自律，是每周堅持刷幾篇最新論文 | PaperDaily #10

相關內容

熱門內容

最新內容

上一篇：量邦科技馮永昌：AI讓量化投資的戰爭升級，交易麵或許不改變其零和博弈的性質｜人工智能研習社

下一篇：深到骨子裏的自律，是每周堅持刷幾篇最新論文 | PaperDaily #10