1. 前言

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

2. 幾個數學概念

1) 梯度（一階導數）

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

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)，這就是所謂的雅可比矩陣：

總結一下,

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 方法中的一種，主要迭代公式如下：

其中，是第 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)

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)

上面例子中由於目標函式是二次凸函式，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)

大概 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, 第一次移動方向為最速下降方向: 
2. 
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)

參考資料：

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