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

時間  2019-11-10
標籤 機器 學習 導數 最優化 方法 基礎 简体版
原文   原文鏈接

 

1. 前言

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

 

2. 幾個數學概念

1) 梯度（一階導數）

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

2) Hesse 矩陣（二階導數）

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

當 f(x) 爲二次函數時，梯度以及 Hesse 矩陣很容易求得。二次函數能夠寫成下列形式：編程

其中 A 是 n 階對稱矩陣，b 是 n 維列向量， c 是常數。f(x) 梯度是 Ax+b, Hesse 矩陣等於 A。數組

3) Jacobi 矩陣

Jacobi 矩陣其實是向量值函數的梯度矩陣，假設F:Rn→Rm 是一個從n維歐氏空間轉換到m維歐氏空間的函數。這個函數由m個實函數組成: 。這些函數的偏導數(若是存在)能夠組成一個m行n列的矩陣(m by n)，這就是所謂的雅可比矩陣：app

總結一下,機器學習

a) 若是 f(x) 是一個標量函數，那麼雅克比矩陣是一個向量，等於 f(x) 的梯度， Hesse 矩陣是一個二維矩陣。若是 f(x) 是一個向量值函數，那麼Jacobi 矩陣是一個二維矩陣，Hesse 矩陣是一個三維矩陣。ide

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

 

3. 優化方法

1) Gradient Descent

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

其中，是第 k 次迭代咱們選擇移動的方向，在 steepest descent 中，移動的方向設定爲梯度的負方向， 是第 k 次迭代用 line search 方法選擇移動的距離，每次移動的距離係數能夠相同，也能夠不一樣，有時候咱們也叫學習率（learning rate）。在數學上，移動的距離能夠經過 line search 令導數爲零找到該方向上的最小值，可是在實際編程的過程當中，這樣計算的代價太大，咱們通常能夠將它設定位一個常量。考慮一個包含三個變量的函數 ，計算梯度獲得 。設定 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)
View Code

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

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

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

2) Newton's method

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

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

第二步，把 x 看作自變量， 全部帶有 x^k 的項看作常量，令一階導數爲 0 ，便可求近似函數的最小值：

即：

第三步，將當前的最小值設定近似函數的最小值（或者乘以步長)。

與 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)
View Code

上面例子中因爲目標函數是二次凸函數，Taylor 展開等於原函數，因此能一次就求出最優解。

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

  1. Hesse 矩陣不可逆時沒法計算
  2. 矩陣的逆計算複雜爲 n 的立方，當問題規模比較大時，計算量很大，解決的辦法是採用擬牛頓法如 BFGS, L-BFGS, DFP, Broyden's Algorithm 進行近似。
  3. 若是初始值離局部極小值太遠，Taylor 展開並不能對原函數進行良好的近似

3) Levenberg–Marquardt Algorithm

Levenberg–Marquardt algorithm 能結合以上兩種優化方法的優勢，並對二者的不足作出改進。與 line search 的方法不一樣，LMA 屬於一種「信賴域法」(trust region)，牛頓法實際上也能夠看作一種信賴域法，即利用局部信息對函數進行建模近似，求取局部最小值。所謂的信賴域法，就是從初始點開始，先假設一個能夠信賴的最大位移 s（牛頓法裏面 s 爲無窮大），而後在以當前點爲中心，以 s 爲半徑的區域內，經過尋找目標函數的一個近似函數（二次的）的最優勢，來求解獲得真正的位移。在獲得了位移以後，再計算目標函數值，若是其使目標函數值的降低知足了必定條件，那麼就說明這個位移是可靠的，則繼續按此規則迭代計算下去；若是其不能使目標函數值的降低知足必定的條件，則應減少信賴域的範圍，再從新求解。

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

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

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

位移是多少時獲得 S 函數的最小值呢？經過幾何的概念，當殘差 垂直於 J 矩陣的 span 空間時， S 取得最小（至於爲何？請參考以前博客的最後一部分)

咱們將這個公式略加修改，加入阻尼係數獲得：

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

算法過程以下：

  1. 給定一個初識值 x0
  2. 當 而且沒有到達最大迭代次數時
  3. 重複執行:
    • 算出移動向量
    • 計算更新值：
    • 計算目標函數真實減小量與預測減小量的比率
    • if ，接受更新值
    • else if ，說明近似效果很好，接受更新值，擴大可信域（即減少阻尼係數）
    • else: 目標函數在變大，拒絕更新值，減少可信域（即增長阻尼係數）
  4. 直到達到最大迭代次數

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

展現了 zig-zagging 鋸齒現象：

用 LMA 優化效率如何。套用到咱們以前 LMA 公式中，有：

代碼以下：

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)
View Code

大概 5 次迭代就能夠獲得最優解 (1, 1).

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

4) Conjugate Gradients

共軛梯度法也是優化模型常常常常要用到的一個方法，背後的數學公式和原理稍微複雜一些，光這一個優化方法就能夠寫一篇很長的博文了，因此這裏並不打算詳細講解每一步的推導過程，只簡單寫一下算法的實現過程。與最速梯度降低的不一樣，共軛梯度的優勢主要體如今選擇搜索方向上。在瞭解共軛梯度法以前，咱們首先簡單瞭解一下共軛方向：

共軛方向和馬氏距離的定義有相似之處，他們都考慮了全局的數據分佈。如上圖，d(1) 方向與二次函數的等值線相切，d(1) 的共軛方向 d(2) 則指向橢圓的中心。因此對於二維的二次函數，若是在兩個共軛方向上進行一維搜索，通過兩次迭代必然達到最小點。前面咱們說過，等值線橢圓的形狀由 Hesse 矩陣決定，那麼，上圖的兩個方向關於 Hessen 矩陣正交，共軛方向的定義以下：

若是橢圓是一個正圓， Hessen 矩陣是一個單位矩陣，上面等價於歐幾里得空間中的正交。

在優化過程當中，若是咱們肯定了移動方向（GD：垂直於等值線，CG：共軛方向），而後在該方向上搜索極小值點（剛好與該處的等值線相切），而後移動到最小值點，重複以上過程，那麼 Gradient Descent 和 Conjugate gradient descent 的優化過程能夠用下圖的綠線與紅線表示：

講了這麼多，共軛梯度算法到底是如何算的呢？

  1. 給定一個出發點 x0 和一箇中止參數 e, 第一次移動方向爲最速降低方向:
  3. while :
    • 用 Newton-Raphson 迭代計算移動距離，以便在該搜索方向移動到極小，公式就不寫了，具體思路就是利用一階梯度的信息向極小值點跳躍搜索
    • 移動當前的優化解 x：
    • 用 Gram-Schmidt 方法構造下一個共軛方向，即 , 按照 的肯定公式又能夠分爲 FR 方法和 PR 和 HS 等。

在不少的資料中，介紹共軛梯度法都舉了一個求線性方程組 Ax = b 近似解的例子，實際上就至關於這裏所說的

仍是用最開始的目標函數     來編寫共軛梯度法的優化代碼：

 

# 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)
View Code

 

  

 

 

參考資料：

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

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程