機器學習中導數最優化方法(基礎篇)

jopen 9年前發布 | 50K 次閱讀機器學習

作者：daniel-D

1. 前言

熟悉機器學習的童鞋都知道，優化方法是其中一個非常重要的話題，最常見的情形就是利用目標函數的導數通過多次迭代來求解無約束最優化問題。實現簡單，coding 方便，是訓練模型的必備利器之一。這篇博客主要總結一下使用導數的最優化方法的幾個基本方法，梳理梳理相關的數學知識，本人也是一邊寫一邊學，如有問題，歡迎指正，共同學習，一起進步。

2. 幾個數學概念

1) 梯度（一階導數）

考慮一座在 (x1, x2) 點高度是 f(x1, x2) 的山。那么，某一點的梯度方向是在該點坡度最陡的方向，而梯度的大小告訴我們坡度到底有多陡。注意，梯度也可以告訴我們不在最快變化方向的其他方向的變化速度（二維情況下，按照梯度方向傾斜的圓在平面上投影成一個橢圓）。對于一個含有 n 個變量的標量函數，即函數輸入一個 n 維的向量，輸出一個數值，梯度可以定義為：

2) Hesse 矩陣（二階導數）

Hesse 矩陣常被應用于牛頓法解決的大規模優化問題(后面會介紹)，主要形式如下：

當 f(x) 為二次函數時，梯度以及 Hesse 矩陣很容易求得。二次函數可以寫成下列形式： $f(x)=\frac{1}{2}x^TAx + b^Tx+c$ 其中 A 是 n 階對稱矩陣，b 是 n 維列向量， c 是常數。f(x) 梯度是 Ax+b, Hesse 矩陣等于 A。

3) Jacobi 矩陣

Jacobi 矩陣實際上是向量值函數的梯度矩陣，假設F:Rn→Rm 是一個從n維歐氏空間轉換到m維歐氏空間的函數。這個函數由m個實函數組成: $y(x) = [y_1(x_1, ..., x_n),\cdots , y_m(x_1, ...,x_n)]$ 。這些函數的偏導數(如果存在)可以組成一個m行n列的矩陣(m by n)，這就是所謂的雅可比矩陣：

總結一下,

a) 如果 f(x) 是一個標量函數，那么雅克比矩陣是一個向量，等于 f(x) 的梯度， Hesse 矩陣是一個二維矩陣。如果 f(x) 是一個向量值函數，那么Jacobi 矩陣是一個二維矩陣，Hesse 矩陣是一個三維矩陣。

b) 梯度是 Jacobian 矩陣的特例，梯度的 jacobian 矩陣就是 Hesse 矩陣（一階偏導與二階偏導的關系）。

3. 優化方法

1) Gradient Descent

Gradient descent 又叫 steepest descent，是利用一階的梯度信息找到函數局部最優解的一種方法，也是機器學習里面最簡單最常用的一種優化方法。Gradient descent 是 line search 方法中的一種，主要迭代公式如下：

$x_{k+1} = x_k + \alpha _k\mathbf{p}_k$ 其中， $\mathbf{p}_k$ 是第 k 次迭代我們選擇移動的方向，在 steepest descent 中，移動的方向設定為梯度的負方向， $\alpha _k$ 是第 k 次迭代用 line search 方法選擇移動的距離，每次移動的距離系數可以相同，也可以不同，有時候我們也叫學習率（learning rate）。在數學上，移動的距離可以通過 line search 令導數為零找到該方向上的最小值，但是在實際編程的過程中，這樣計算的代價太大，我們一般可以將它設定位一個常量。考慮一個包含三個變量的函數 $f(\mathbf{x})=0.5x_1^2+0.2x_2^2+0.6x_3^2$ ，計算梯度得到 $\nabla f(\mathbf{x})=(x_1,0.4x_2,1.2x_3)$ 。設定 learning rate = 1，算法代碼如下：

# Code from Chapter 11 of Machine Learning: An Algorithmic Perspective
# by Stephen Marsland (http://seat.massey.ac.nz/personal/s.r.marsland/MLBook.html)

# Gradient Descent using steepest descent

from numpy import *

def Jacobian(x):
    return array([x[0], 0.4*x[1], 1.2*x[2]])

def steepest(x0):

    i = 0 
    iMax = 10
    x = x0
    Delta = 1
    alpha = 1

    while i<iMax and Delta>10**(-5):
        p = -Jacobian(x)
        xOld = x
        x = x + alpha*p
        Delta = sum((x-xOld)**2)
        print 'epoch', i, ':'
        print x, '\n'
        i += 1

x0 = array([-2,2,-2])
steepest(x0)

Steepest gradient 方法得到的是局部最優解，如果目標函數是一個凸優化問題，那么局部最優解就是全局最優解，理想的優化效果如下圖，值得注意一點的是，每一次迭代的移動方向都與出發點的等高線垂直：

需要指出的是，在某些情況下，最速下降法存在 鋸齒現象（ zig-zagging）將會導致收斂速度變慢:

粗略來講，在二次函數中，橢球面的形狀受 hesse 矩陣的條件數影響，長軸與短軸對應矩陣的最小特征值和最大特征值的方向，其大小與特征值的平方根成反比，最大特征值與最小特征值相差越大，橢球面越扁，那么優化路徑需要走很大的彎路，計算效率很低。

2) Newton’s method

在最速下降法中，我們看到，該方法主要利用的是目標函數的局部性質，具有一定的“盲目性”。牛頓法則是利用局部的一階和二階偏導信息，推測整個目標函數的形狀，進而可以求得出近似函數的全局最小值，然后將當前的最小值設定近似函數的最小值。相比最速下降法，牛頓法帶有一定對全局的預測性，收斂性質也更優良。牛頓法的主要推導過程如下：

第一步，利用 Taylor 級數求得原目標函數的二階近似：

$f(\mathbf{x})\approx \phi (\mathbf{x})=f(\mathbf{x}^{(k)})+\nabla f(\mathbf{x}^{(k)})(\mathbf{x}-\mathbf{x}^{(k)})+\frac{1}{2}(\mathbf{x}-\mathbf{x}^{(k)})^T\nabla ^2f(\mathbf{x}^{(k)})(\mathbf{x}-\mathbf{x}^{(k)})$ 第二步，把 x 看做自變量，所有帶有 x^k 的項看做常量，令一階導數為 0 ，即可求近似函數的最小值： $\mathbf{p}_k = -(\nabla ^2f(\mathbf{x_k}))^{-1}\nabla f(\mathbf{x_k})$ 即： $\mathbf{p}_k = -Hesse^{-1}Jacobi$ 第三步，將當前的最小值設定近似函數的最小值（或者乘以步長)。

與 1) 中優化問題相同，牛頓法的代碼如下：

Newton.py

# Code from Chapter 11 of Machine Learning: An Algorithmic Perspective
# by Stephen Marsland (http://seat.massey.ac.nz/personal/s.r.marsland/MLBook.html)

# Gradient Descent using Newton's method
from numpy import *

def Jacobian(x):
    return array([x[0], 0.4*x[1], 1.2*x[2]])

def Hessian(x):
    return array([[1,0,0],[0,0.4,0],[0,0,1.2]])

def Newton(x0):

    i = 0
    iMax = 10
    x = x0
    Delta = 1
    alpha = 1

    while i<iMax and Delta>10**(-5):
        p = -dot(linalg.inv(Hessian(x)),Jacobian(x))
        xOld = x
        x = x + alpha*p
        Delta = sum((x-xOld)**2)
        i += 1
    print x

x0 = array([-2,2,-2])
Newton(x0)

上面例子中由于目標函數是二次凸函數，Taylor 展開等于原函數，所以能一次就求出最優解。

牛頓法主要存在的問題是：

Hesse 矩陣不可逆時無法計算
矩陣的逆計算復雜為 n 的立方，當問題規模比較大時，計算量很大，解決的辦法是采用擬牛頓法如 BFGS, L-BFGS, DFP, Broyden’s Algorithm 進行近似。
如果初始值離局部極小值太遠，Taylor 展開并不能對原函數進行良好的近似

3) Levenberg–Marquardt Algorithm

Levenberg–Marquardt algorithm 能結合以上兩種優化方法的優點，并對兩者的不足做出改進。與 line search 的方法不同，LMA 屬于一種“信賴域法”(trust region)，牛頓法實際上也可以看做一種信賴域法，即利用局部信息對函數進行建模近似，求取局部最小值。所謂的信賴域法，就是從初始點開始，先假設一個可以信賴的最大位移 s（牛頓法里面 s 為無窮大），然后在以當前點為中心，以 s 為半徑的區域內，通過尋找目標函數的一個近似函數（二次的）的最優點，來求解得到真正的位移。在得到了位移之后，再計算目標函數值，如果其使目標函數值的下降滿足了一定條件，那么就說明這個位移是可靠的，則繼續按此規則迭代計算下去；如果其不能使目標函數值的下降滿足一定的條件，則應減小信賴域的范圍，再重新求解。

LMA 最早提出是用來解決最小二乘法曲線擬合的優化問題的，對于隨機初始化的已知參數 beta，求得的目標值為：

對擬合曲線函數進行一階 Jacobi 矩陣的近似：

進而推測出 S 函數的周邊信息：

位移是多少時得到 S 函數的最小值呢？通過幾何的概念，當殘差 $\mathbf{y}-f(\beta)-J\delta$ 垂直于 J 矩陣的 span 空間時， S 取得最小（至于為什么？請參考之前博客的最后一部分)

我們將這個公式略加修改，加入阻尼系數得到：

就是萊文貝格－馬夸特方法。這種方法只計算了一階偏導，而且不是目標函數的 Jacobia 矩陣，而是擬合函數的 Jacobia 矩陣。當 $\lambda$ 大的時候可信域小，這種算法會接近最速下降法， $\lambda$ 小的時候可信域大，會接近高斯-牛頓方法。

算法過程如下：

給定一個初識值 x0
當 $\mathbf{J}^T[\mathbf{y-f(\beta)}] > tolerance$ 并且沒有到達最大迭代次數時
重復執行:
- 算出移動向量 $\delta$
- 計算更新值： $\mathbf{x_{new}}=\mathbf{x}+\delta$
- 計算目標函數真實減少量與預測減少量的比率 $\rho$
- if $0 < \rho < 0.25$ ，接受更新值
- else if $\rho > 0.25$ ，說明近似效果很好，接受更新值，擴大可信域（即減小阻尼系數）
- else: 目標函數在變大，拒絕更新值，減小可信域（即增加阻尼系數）
直到達到最大迭代次數

維基百科在介紹 Gradient descent 時用包含了細長峽谷的 Rosenbrock function

展示了 zig-zagging 鋸齒現象：

用 LMA 優化效率如何。套用到我們之前 LMA 公式中，有： ${ y }\quad =\quad 0\\ f(\beta )\quad =\quad (10(x_{ 2 }-x_{ 1 }^{ 2 }),\quad 1-x_{ 1 })^{ T }\\ J\quad =\quad \begin{pmatrix} -20{ x }_{ 1 } & 10 \\ -1 & 0 \end{pmatrix}$ 代碼如下：

LevenbergMarquardt.py

# Code from Chapter 11 of Machine Learning: An Algorithmic Perspective
# by Stephen Marsland (http://seat.massey.ac.nz/personal/s.r.marsland/MLBook.html)

# The Levenberg Marquardt algorithm
from numpy import *

def function(p):
    r = array([10*(p[1]-p[0]**2),(1-p[0])])
    fp = dot(transpose(r),r) #= 100*(p[1]-p[0]**2)**2 + (1-p[0])**2
    J = (array([[-20*p[0],10],[-1,0]]))
    grad = dot(transpose(J),transpose(r))
    return fp,r,grad,J

def lm(p0,tol=10**(-5),maxits=100):

    nvars=shape(p0)[0]
    nu=0.01
    p = p0
    fp,r,grad,J = function(p)
    e = sum(dot(transpose(r),r))
    nits = 0
    while nits<maxits and linalg.norm(grad)>tol:
        nits += 1
        fp,r,grad,J = function(p)
        H=dot(transpose(J),J) + nu*eye(nvars)

        pnew = zeros(shape(p))
        nits2 = 0
        while (p!=pnew).all() and nits2<maxits:
            nits2 += 1
            dp,resid,rank,s = linalg.lstsq(H,grad)
            pnew = p - dp
            fpnew,rnew,gradnew,Jnew = function(pnew)
            enew = sum(dot(transpose(rnew),rnew))
            rho = linalg.norm(dot(transpose(r),r)-dot(transpose(rnew),rnew))
            rho /= linalg.norm(dot(transpose(grad),pnew-p))

            if rho>0:
                update = 1
                p = pnew
                e = enew
                if rho>0.25:
                    nu=nu/10
            else: 
                nu=nu*10
                update = 0
        print fp, p, e, linalg.norm(grad), nu

p0 = array([-1.92,2])
lm(p0)

大概 5 次迭代就可以得到最優解 (1, 1).

Levenberg–Marquardt algorithm 對局部極小值很敏感，維基百科舉了一個二乘法曲線擬合的例子，當使用不同的初始值時，得到的結果差距很大，我這里也有 python 代碼，就不細說了。

4) Conjugate Gradients

共軛梯度法也是優化模型經常經常要用到的一個方法，背后的數學公式和原理稍微復雜一些，光這一個優化方法就可以寫一篇很長的博文了，所以這里并不打算詳細講解每一步的推導過程，只簡單寫一下算法的實現過程。與最速梯度下降的不同，共軛梯度的優點主要體現在選擇搜索方向上。在了解共軛梯度法之前，我們首先簡單了解一下共軛方向：

共軛方向和馬氏距離的定義有類似之處，他們都考慮了全局的數據分布。如上圖，d(1) 方向與二次函數的等值線相切，d(1) 的共軛方向 d(2) 則指向橢圓的中心。所以對于二維的二次函數，如果在兩個共軛方向上進行一維搜索，經過兩次迭代必然達到最小點。前面我們說過，等值線橢圓的形狀由 Hesse 矩陣決定，那么，上圖的兩個方向關于 Hessen 矩陣正交，共軛方向的定義如下： $\mathbf{d}^{(1)T}\mathbf{Ad}^{(2)} = 0}$ 如果橢圓是一個正圓， Hessen 矩陣是一個單位矩陣，上面等價于歐幾里得空間中的正交。

在優化過程中，如果我們確定了移動方向（GD：垂直于等值線，CG：共軛方向），然后在該方向上搜索極小值點（恰好與該處的等值線相切），然后移動到最小值點，重復以上過程，那么 Gradient Descent 和 Conjugate gradient descent 的優化過程可以用下圖的綠線與紅線表示：

講了這么多，共軛梯度算法究竟是如何算的呢？

給定一個出發點 x0 和一個停止參數 e, 第一次移動方向為最速下降方向: $\mathbf{p}_0 = -\nabla f(\mathbf{x})}$
$\mathbf{p}_{new} = \mathbf{p}_0$
while :
- 用 Newton-Raphson 迭代計算移動距離，以便在該搜索方向移動到極小，公式就不寫了，具體思路就是利用一階梯度的信息向極小值點跳躍搜索
- 移動當前的優化解 x： $\mathbf{x}_{new} = \mathbf{x} + \alpha _k\mathbf{p}$
- 用 Gram-Schmidt 方法構造下一個共軛方向，即 $\mathbf{p}_{new} = \nabla f(\mathbf{x}_{new}) + \beta _{k+1}\mathbf{p}$ , 按照 $\beta$ 的確定公式又可以分為 FR 方法和 PR 和 HS 等。

在很多的資料中，介紹共軛梯度法都舉了一個求線性方程組 Ax = b 近似解的例子，實際上就相當于這里所說的 $\min f(\mathbf{x}) = \frac{1}{2}||\mathbf{Ax-b}||^2$

還是用最開始的目標函數 $f(\mathbf{x})=0.5x_1^2 + 0.2x_2^2 + 0.6x_3^2$ 來編寫共軛梯度法的優化代碼：

# Code from Chapter 11 of Machine Learning: An Algorithmic Perspective
# by Stephen Marsland (http://seat.massey.ac.nz/personal/s.r.marsland/MLBook.html)

# The conjugate gradients algorithm

from numpy import *

def Jacobian(x):
    #return array([.4*x[0],2*x[1]])
    return array([x[0], 0.4*x[1], 1.2*x[2]])

def Hessian(x):
    #return array([[.2,0],[0,1]])
    return array([[1,0,0],[0,0.4,0],[0,0,1.2]])

def CG(x0):

    i=0
    k=0

    r = -Jacobian(x0)
    p=r

    betaTop = dot(r.transpose(),r)
    beta0 = betaTop

    iMax = 3
    epsilon = 10**(-2)
    jMax = 5

    # Restart every nDim iterations
    nRestart = shape(x0)[0]
    x = x0

    while i < iMax and betaTop > epsilon**2*beta0:
        j=0
        dp = dot(p.transpose(),p)
        alpha = (epsilon+1)**2
        # Newton-Raphson iteration
        while j < jMax and alpha**2 * dp > epsilon**2:
            # Line search
            alpha = -dot(Jacobian(x).transpose(),p) / (dot(p.transpose(),dot(Hessian(x),p)))
            print "N-R",x, alpha, p
            x = x + alpha * p
            j += 1
        print x
        # Now construct beta
        r = -Jacobian(x)
        print "r: ", r
        betaBottom = betaTop
        betaTop = dot(r.transpose(),r)
        beta = betaTop/betaBottom
        print "Beta: ",beta
        # Update the estimate
        p = r + beta*p
        print "p: ",p
        print "----"
        k += 1

        if k==nRestart or dot(r.transpose(),p) <= 0:
            p = r
            k = 0
            print "Restarting"
        i +=1

    print x

x0 = array([-2,2,-2])
CG(x0)

參考資料：

[1] Machine Learning: An Algorithmic Perspective, chapter 11
[2] 最優化理論與算法（第2版），陳寶林
[3] wikipedia

文章出處：http://www.cnblogs.com/daniel-D/

本文由用戶 jopen 自行上傳分享，僅供網友學習交流。所有權歸原作者，若您的權利被侵害，請聯系管理員。

轉載本站原創文章，請注明出處，并保留原始鏈接、圖片水印。

本站是一個以用戶分享為主的開源技術平臺，歡迎各類分享！

本文地址：http://www.baiduhome.net/lib/view/open1425889892388.html

機器學習