Optimization of Machine Learning

機器學習就是需要找到模型的鞍點，也就是最優點。因為模型很多時候并不是完全的凸函數，是以如果沒有好的優化方法可能會跑不到極值點，或者是局部極值，甚至是偏離。是以選擇一個良好的優化方法是至關重要的。首先是比較正常的優化方法：梯度下降。以下介紹的這些算法都不是用于當個算法，可以試用于能可微的所有算法。

Gradient Descent

常見會用在logistics regression或者是linear regression裡面。比如logistics regression，這個模型直接等于0是求不出解析解的，所有隻能用一些疊代方法求解。當樣本的label是

每一個樣本正确機率：

最大似然函數

化簡：

target function就是如上的

，要做的就是優化上面的損失函數。

這個就是目前點的梯度，而函數的更新方向就是向着梯度的負方向進行的，是以最後是要減：

最優點很明顯是在左邊，而待優化點是

點，而在

點的梯度是上升的，是以負方向才是正确的，最後得到的梯度就是要減去了。這個就是gradient descent。對于步長

，可以選擇

，t就是循環次數，随着疊代次數的增多，越來越靠近極值點，再相同的疊代步長是有可能越過極值點的。事實上步長的選擇也是一個很重要的方法，後面會使用到。梯度下降是用給一個平面去拟合，也就是一次曲線去拟合，因為梯度下降是用泰勒一階展開的，是以是一次曲線。這樣造成的問題就是，gradient descent方法的視野很小，因為它隻有一階拟合，是以目前找到那個方向并不一定是全局最優的方向，可能是有點偏離，是以經常出現之字形的走勢，疊代的速度也不快。它隻會看目前一階導數哪個最适合，而不會去管二次或者三次的。

梯度下降有三種：全量梯度下降，所有是資料一起做疊代更新。資料量很大的話速度回很慢。小批量梯度下降，一部分一部分資料的就行疊代資料。随機梯度下降SGD，随機選取幾個進行疊代，可能疊代的方向會有偏差，但是随着時間流逝大方向還是一樣的。代碼實作前面的logistics regression中已經有了，不再重複。

給予這種不足，給出了牛頓法的改進。

NewTon Method

既然梯度下降做的是一階拟合，那麼試一下二階拟合。泰勒展開二階：

兩邊對

求導數：

代換一下：

牛頓法有點像坐标上升的方法，先疊代一個坐标，再疊代另一個，和SMO類似。可以看到牛頓法這個時候就是二次曲線的拟合了：

左邊的就是成功拟合的，右邊的那個就是方向錯誤，拟合沒有成功的例子。是以牛頓法不一定會按照正确的方向拟合。上面牛頓法的式子是對于單變量的，如果是對變量，那麼下面的二階導要用到Hession矩陣。是以對于多變量的牛頓法：

稱為是牛頓方向，剛剛提到，疊代方向要和一階導數方向相反，也就是

，這個想法和剛剛提出的沿着梯度方向的反方向下降的思路是一樣的。分解一下這個公式：

，是以Hession矩陣要是正定矩陣而且非奇異，而正定矩陣這個條件是很強的，不是在特定條件下達不到，是以牛頓法雖然下降的速度可以很快，但是方向不一定是正确的，是以牛頓法要使用的話一定要要求是靠近極值點的情況下使用，否則容易偏離。從本質上去看，牛頓法是二階收斂，梯度下降是一階收斂，是以牛頓法就更快。如果更通俗地說的話，比如你想找一條最短的路徑走到一個盆地的最底部，梯度下降法每次隻從你目前所處位置選一個坡度最大的方向走一步，牛頓法在選擇方向時，不僅會考慮坡度是否夠大，還會考慮你走了一步之後，坡度是否會變得更大。是以，可以說牛頓法比梯度下降法看得更遠一點，能更快地走到最底部。（牛頓法目光更加長遠，是以少走彎路；相對而言，梯度下降法隻考慮了局部的最優，沒有全局思想。）這一段簡單來說就是梯度下降僅僅考慮了梯度，而牛頓法不僅考慮了梯度，而且還考慮了梯度的變化，也就是梯度的梯度。

代碼實作

def sigmoid(
        x1, x2, theta1, theta2, theta3):
    z = (theta1*x1+ theta2*x2+theta3).astype("float_")
    return 1.0 / (1.0 + np.exp(-z))

'''cross entropy loss function
'''
def cost(
        x1, x2, y, theta_1, theta_2, theta_3):
    sigmoid_probs = sigmoid(x1,x2, theta_1, theta_2,theta_3)
    return -np.mean(y * np.log(sigmoid_probs)
                  + (1 - y) * np.log(1 - sigmoid_probs))
           

sigmoid函數和cost函數，使用牛頓法代替梯度下降做邏輯回歸。

def gradient(
        x_1, x_2, y, theta_1, theta_2, theta_3):
    sigmoid_probs = sigmoid(x_1,x_2, theta_1, theta_2,theta_3)
    return np.array([np.sum((y - sigmoid_probs) * x_1),
                     np.sum((y - sigmoid_probs) * x_2),
                     np.sum(y - sigmoid_probs)])

           

一階導數，資料就是二維的，再加上一個偏置值就是三個導數了。

'''calculate hassion matrix'''
def Hassion(
        x1, x2, y, theta1, theta2, theta3):
    sigmoid_probs = sigmoid(x1,x2, theta1,theta2,theta3)
    d1 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x1 * x1)
    d2 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x1 * x2)
    d3 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x1 )
    d4 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x2 * x1)
    d5 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x2 * x2)
    d6 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x2 )
    d7 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x1 )
    d8 = np.sum((sigmoid_probs * (1 - sigmoid_probs)) * x2 )
    d9 = np.sum((sigmoid_probs * (1 - sigmoid_probs)))
    H = np.array([[d1, d2,d3],[d4, d5,d6],[d7,d8,d9]])
    return H
           

方法比較笨，直接一個一個疊代。

'''update parameter by newton method'''
def newton_method(x1, x2, y):
    #initialize parameter
    theta1 = 0.001
    theta2 = -0.4
    theta3 = 0.6
    delta_loss = np.Infinity
    loss = cost(x1, x2, y, theta1, theta2, theta3)
    kxi = 0.0000001
    max_iteration = 10000
    i = 0
    print('while i = ', i, ' Loss = ', loss, 'delta_loss = ', delta_loss)
    while np.abs(delta_loss) > kxi and i < max_iteration:
        i += 1
        g = gradient(x1, x2, y, theta1, theta2, theta3)
        hassion = Hassion(x1, x2, y, theta1, theta2, theta3)
        H_inverse = np.linalg.inv(hassion)
        delta = H_inverse @ g.T
        delta_theta_1 = delta[0]
        delta_theta_2 = delta[1]
        delta_theta_3 = delta[2]
        theta1 += delta_theta_1
        theta2 += delta_theta_2
        theta3 += delta_theta_3
        loss_new = cost(x1, x2, y, theta1, theta2, theta3)
        delta_loss = loss - loss_new
        loss = loss_new
        print('while i = ', i, ' Loss = ', loss, 'delta_loss = ', delta_loss)
    return np.array([theta1, theta2, theta3])
           

運作函數。

def get_param_target(self):
        np.random.seed(12)
        num_observations = 5000
        x1 = np.random.multivariate_normal([0, 0], [[1, .75], [.75, 1]], num_observations)
        x2 = np.random.multivariate_normal([1, 4], [[1, .75], [.75, 1]], num_observations)
        simulated_separableish_features = np.vstack((x1, x2)).astype(np.float32)
        simulated_labels = np.hstack((np.zeros(num_observations),
                                      np.ones(num_observations)))
        return (simulated_separableish_features[:, 0], simulated_separableish_features[:, 1], simulated_labels)

           

生成資料。

如果是随機值的話，多試幾次有可能方向是錯誤的。這裡能跑這麼快是因為初始值設定的好。

對于牛頓法的這個問題還是有改進方法的。對于

步長的算法研究是可以解決這個問題。既然有時候牛頓法下降的方向或許是不對的，可以使用一個阻力來阻擋一下，減慢速度，等到下一次可以重新選擇。是以接下來就是對于

步長的研究。

線搜尋

步長最小化意味着

，

就是方向了,這裡的

都是常數，是以可以看做是一個

的函數。

二分搜尋算法

假設

,如果是一個凸函數，需要

，因為

，p就目前的梯度下降的方向，是以向着梯度的反方向得：

，假如我們知道一個

，那麼就可以知道肯定有一個

使得

，是以二分搜尋的流程很簡單：

①k = 0，

②

，如果

如果

，否則就繼續循環。

二分搜尋也叫精确搜尋，因為這樣不斷的疊代下去是可以很精确的搜尋到一個合适的值，使用二分查找法來求步長的計算複雜度很高, 因為在最小化f(x)f(x)的每次疊代中我們都需要執行一次線搜尋, 而每次線搜尋都要用上述的二分查找算法，因為這個精确的值在一些凸函數裡面會是很小很小的，甚至後面多加一個小數，這樣算起來複雜度非常高的。我們可以在犧牲一定的精度的條件下來加快計算速度, 回溯線搜尋是一種近似線搜尋算法。

回溯線搜尋

對于上訴二分法的問題，這裡提出了改進，犧牲掉一些精度，求一個大概的範圍就好，于是有：

上訴條件稱為充分下降條件，

這兩個準則就是Armijo-Goldstein準則，①目标函數值應該有足夠的下降；②一維搜尋的步長α不應該太小。這兩個目标其實很明顯，足夠的下降其實就是式子1：

可以看到接受這個條件有兩個區間，有時候會選擇到第一個區間的内容，也就是第一個區間的内容，是以第二條式子的作用就是舍得步長不要太小了。

為什麼

①對于

的選擇是一定要在(0,0.5)之間的，如果沒有這個條件，會影響算法的超線性收斂性。這個怎麼收斂的我就不知道了，還得看看其他的凸優化書籍。

做泰勒展開，

，去掉高階無窮小，

，和上式就是相差了一個

而已，這樣就保證了小于的這個條件。

③由于

，是以

，而且

，是以有

，這就保證了a不能太小

綜上，這就得到了Armijo-Goldstein準則。後面用的Armijo搜尋隻是用了第一條式子進行搜尋。

阻尼牛頓法

回到剛剛的牛頓法，提到有可能會存在Hession Matrix不是正定的情況，我們可以加一點阻礙，也就是加上一個步長的限制，這個步長就是由Arimji搜尋得到，這裡隻是會用到第一條準則，第二三條先不用。

class DampedNewton(object):

    def __init__(self, feature, label, iterMax, sigma, delta):
        self.feature = feature
        self.label = label
        self.iterMax = iterMax
        self.sigma = sigma
        self.delta = delta
        self.w = None

    def get_error(self, w):
        return (self.label - self.feature*w).T * (self.label - self.feature*w)/2

    def first_derivative(self):
        m, n = np.shape(self.feature)
        g = np.mat(np.zeros((n, 1)))
        for i in range(m):
            err = self.label[i, 0] - self.feature[i, ]*self.w
            for j in range(n):
                g[j, ] -= err*self.feature[i, j]
        return g

    def second_derivative(self):
        m, n = np.shape(self.feature)
        G = np.mat(np.zeros((n, n)))
        for i in range(m):
            x_left = self.feature[i, ].T
            x_right = self.feature[i, ]
            G += x_left * x_right
        return G

    def get_min_m(self, d, g):
        m = 0
        while True:
            w_new = self.w + pow(self.sigma, m)*d
            left = self.get_error(w_new)
            right = self.get_error(self.w) + self.delta*pow(self.sigma, m)*g.T*d
            if left <= right:
                break
            else:
                m += 1
        return m

    def newton(self):
        n = np.shape(self.feature)[1]
        self.w = np.mat(np.zeros((n, 1)))
        it = 0
        while it <= self.iterMax:
            g = self.first_derivative()
            G = self.second_derivative()
            d = -G.I * g
            m = self.get_min_m(d, g)
            self.w += pow(self.sigma, m) * d
            if it % 100 == 0:
                print('it: ', it, ' error: ', self.get_error(self.w)[0, 0])
            it += 1
        return self.w

           

求海塞矩陣換了一個寫法，好看一點，但是效果都是一樣的。

線性回歸的損失函數：

線性回歸的一階導數：

線性回歸的二階導數：

對照上面的公式是可以一一對應的。對于Armrji搜尋隻是用了第一條公式：

一旦超過這個限制就可以退出了。

效果其實還是可以的。當然，也可以三條準則都使用上，隻不過使用一條一般也足夠了。

到這裡線搜尋還沒有完。

Wolfe-Powell準則

有Armijo-Goldstein準則還是不夠的，如下圖：

可以看到(a,b),(c,d)就是選擇的區間，但是很明顯這些區間已經避開了最低的點，當然這不是一定會，但是有這個可能，為了解決這個問題就出現了Wolfe-Powell準則。

Wolfe-Powell準則也有兩個表達式，第一個表達式就是

第二條式子就去掉了，因為它會很容易的把極值點排除在外，是以并不需要這個條件。添加的就是：

這條式子的直覺解釋就是，目前可選點的斜率要大于初始點斜率的

倍。，而初始點（ α=0 的點）處的切線是比 e 點處的切線要“斜”的，由于

∈(ρ,1) ，使得 e 點處的切線變得“不那麼斜”了：

e和f後面的點就是滿足的了，

這個條件是沒有的了。這個條件還是不太好，有時候會給出更嚴謹的一個條件：

，其實就是兩邊都給絕對值，但是右邊的本來就是一個負數了，加上符号就是正數。這個條件就更強了：

直接限定在了(e,g)和(f,h)範圍之内了。這就是整個Wolfe-Powell準則。

線搜尋的總結

線搜尋到這裡基本就結束了。尋找一個最佳步長就是要求優化之後的函數值是要最小的，于是給出的目标函數就是：

。解決這個問題有兩種方法，一種就是二分法的查找，這種方法可能非常的精确搜尋到一個最佳的值，但是他的計算複雜度有點高；有時候我們其實不太需要太過精确的值，我們隻是需要一個大概模糊的就好了，于是出現了回溯搜尋，這是屬于模糊搜尋，給出的就是Armijo-Goldstein準則，但是這兩個準則也有不好的地方，有時候會把極值點排除在外，于是引入曲率條件，把Armijo-Goldstein準則的第二條式子換掉，就出現了Wolfe-Powell準則，條件再強一點兩邊加上絕對值，這樣就出現了最終的Wolfe-Powell準則。

對于步長的研究并不針對某一個算法，對于優化下降的算法都可以使用，梯度下降，牛頓法，拟牛頓法都可以用到，具有很強的普适性。梯度下降的步長也是可以通過這種方式進行選擇最優步長，牛頓法用Armijo搜尋的方法是可以得到全局牛頓法，也叫阻尼牛頓法，這樣可以使得疊代方向可以避免向錯誤的方向進行，增加點阻力。

這篇文章主要的研究對象還是牛頓法，是以下面的三個算法分别就是DFP，BFGS，LBFGS。

DFP

前面的阻尼牛頓法解決的就是對于疊代方向的問題，但是對于計算複雜度這個問題還沒有得到解決，因為矩陣求拟是一個很大的工程量，如果次元一多是計算複雜度是很大的，是以拟牛頓法基本上都是構造一個和Hession矩陣相似的矩陣來替代。但是問題來了，用近似矩陣替代，效果可能會不好，給出如下證明：

首先由牛頓法可以得到：

由于

，是以右邊也是要求小于0，那麼自然就是大于0了，這也就證明了Hession矩陣的正定性。在逐漸尋找最優解的過程中要求函數值是下降的，但是有時候Hession矩陣不一定是半正定的，可能會使得函數值不降反升。 拟牛頓法可以使目标函數值沿下降方向走下去，并且到了最後，在極小值點附近，可使構造出來的矩陣與Hesse矩陣“很像”了，這樣，拟牛頓法也會具有牛頓法的二階收斂性。

推導證明

涉及到Hession矩陣，自然就涉及到二階導數矩陣了，Hession矩陣隻是對于多變量二階導數的一種描述方式。Taylor展開：

還是用二次導數做近似，很多函數都可以用二次導數做近似。去掉高次項兩邊求導數：

上面那個方程就是拟牛頓方程。上面那個矩陣就是一個近似矩陣。近似矩陣有很多種，如果按照上面那種方式進行計算，複雜度沒有降多少，是以比較常見的方法就是疊代計算，比較常見的疊代方式：

那麼接下來問題來了，

要這麼設計。一般給出的設計就是

因為希望每次疊代h能有一個微小的變化而不是巨大的變化，這樣才有可能收斂。而且這個結構設計的也很簡單，也是對稱的。

代入上面的拟牛頓方程：

仔細看一下這個式子，右邊的

是一個nx1的向量，而

均為實數，也就是點積。可以直接假設

，于是有

。

代入上面的式子：

一開始我看到這個推導過程，我是一臉懵逼的。原來數學家也有猜的時候，這個過程和前面假設H的形式在凸優化裡面，是用"很顯然"來描述的，嗯，很顯然！我就是他媽的沒看出來顯然在哪了，後面的LBFGS也是，真的很顯然，但是這個我是真的不知道。看到網上的很多解釋都離不開特殊情況，也就是為0或者是直接假設

的，都差不多的解釋。

已知初始正定矩陣H，從一個初始點開始（疊代），用式子

來計算出下一個搜尋方向，并在該方向上求出可使目标函數極小化的步長α，然後用這個步長，将目前點挪到下一個點上，并檢測是否達到了程式中止的條件，如果沒有達到，則用上面所說的[13]式的方法計算出下一個修正矩陣H，并計算下一個搜尋方向……周而複始，直到達到程式中止條件。

是以整一個流程：

①給定一個初值

②搜尋方向：

③利用搜尋得到步長

，可以使用上面提到的Armrji搜尋或者等等的改進方法。

④更新一波

⑤計算

⑥

，轉回去繼續更新。

DFP還是有缺點的：

具體主要内容就是說如果Hession矩陣的是錯誤估計的，你們BFGS會在很少的疊代回到正确的方向，但是DFP在這方面的效果不明顯。至于為什麼不明顯，我從公式還不能直接看出來。既然提到了BFGS，下面就是BFGS了。

BFGS

BFGS和DFP其實是屬于對偶的解法，一個直接求海賽矩陣逆矩陣（DFP），另一個就是求海賽矩陣（BFGS）。還是使用拟牛頓公式：

，而DFP使用的是

，直接求出來的就是逆矩陣了。推導過程其實很簡單，和上面DFP類似，差别不大，甚至除了換一下沒什麼差别的。

假設都和上面DFP的一樣。還是

，E就是矯正矩陣。E的形式和之前的是一樣的

。代入上面的公式：

同樣假設

，同樣選取特殊情況，

，代進矯正函數的表達式：

仔細看一些這個式子，和DFP那個式子是換了一下位置而已，是以這個過程和步驟真的沒有什麼好說的。按照正常套路，一般是這樣：

然而，如果是這樣，複雜度還是存在的，還是得求個導數啊。是以這樣的做法自然是不行的，因為拟牛頓法的改進有一個很重要的誘因就是計算複雜度的問題了。

對于求逆矩陣，有一個非常重要的公式——Sherman-Morrison公式：

用Sherman-Morrison公式改造一下：

至于是怎麼改造的，有興趣自己看吧，原諒我并沒有看懂，懂了我也想不出來。

根據這個結果改造一下：

完美的達到了降低複雜度的目标。但是我還是不知道BFGS為什麼相對DFP來說對于糾正方向要更快一點。這兩個是對偶式子，難到是因為用了Sherman-Morrison公式嗎？

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from DataProcesser import DataProcesser

def bfgs(feature, label, lam, maxCycle):
    n = np.shape(feature)[1]
    w0 = np.mat(np.zeros((n, 1)))
    rho = 0.55
    sigma = 0.4
    Bk = np.eye(n)
    k = 1
    while (k < maxCycle):
        print('Iterator: ', k, ' error: ', get_error(feature, label, w0))
        gk = get_gradient(feature, label, w0, lam)
        dk = np.mat(-np.linalg.solve(Bk, gk))
        m = 0
        mk = 0
        while (m < 20):
            newf = get_result(feature, label, (w0 + rho**m*dk), lam)
            oldf = get_result(feature, label, w0, lam)
            if (newf < oldf + sigma * (rho ** m)*(gk.T*dk)[0, 0]):
                mk = m
                break
            m += 1
        #BFGS
        w = w0 + rho**mk * dk
        sk = w-w0
        yk = get_gradient(feature, label, w, lam) - gk
        if (yk.T * sk > 0):
            Bk = Bk - (Bk * sk * sk.T * Bk) / (sk.T * Bk * sk) + (yk * yk.T) / (yk.T * sk)
        k = k + 1
        w0 = w
    return w0

def get_error(feature, label, w):
    return (label - feature * w).T*(label - feature * w)/2

def get_gradient(feature, label, w, lam):
    err = (label - feature * w).T
    left = err * (-1) * feature
    return left.T + lam * w

def get_result(feature, label, w, lam):
    left = (label - feature * w).T * (label - feature * w)
    right = lam * w.T * w
    return (left + right)/2

           

最主要的部分也就是公式那部分了，沒有什麼好講的。實作效果：

L-BFGS

L-BFGS的出現還是為了節省，但這次是為了節省記憶體。每一次存儲高緯度的資料需要的空間太多了，是以LBFGS的主要作用就是減少記憶體的開銷。不再完整的存儲整一個矩陣了，它對BFGS算法做了近似，而是存儲計算時所需要的

，然後利用哪兩個向量的計算來近似代替。對于

也不是存儲所有的，而是隻儲存那麼後面的m個，也就是最新的m個就好。這樣就把空間複雜度

降到

毫無疑問，LBFGS是從BFGS演化而來的。對于逆矩陣的公式：

，改造一下上式：

假設：

遞推一下：

根據推論，一般的有：

可以看的出，目前的海賽矩陣的逆矩陣是可以由

給出的，是以直接存儲對應的

即可，上面提到過為了壓縮記憶體，可以隻是存儲後面的m個，是以更新一下公式：

公式長是有點長，但是已經有大神給出了一個很好的算法解決這個問題，two-loop-recursion算法：

end for

這個式子到底行不行呢？證明一下理論：

這樣就證明這個算法的正确性。然而其實我根本不關心這個算法正确性，我隻是想知道

這是怎麼想出來的，說實話第一眼看根本沒有get到這個算法就是實作了LBFGS，是以如果有大神知道麻煩私信我！渣渣感激不盡。

按照上面算法得到方向之後就可以使用線搜尋得到合适的步長最後結合更新了。

前面求梯度都一樣的，就是後面的更新有不同：

def lbfgs(feature, label, lam, maxCycle, m = 10):
    n = np.shape(feature)[1]
    w0 = np.mat(np.zeros((n, 1)))
    rho = 0.55
    sigma = 0.4

    H0 = np.eye(n)
    s = []
    y = []

    k = 1
    gk = get_gradient(feature, label, w0, lam)
    dk = -H0 * gk

    while (k < maxCycle):
        print('iter: ', k, ' error:', get_error(feature, label, w0))
        m1 = 0
        mk = 0
        gk = get_gradient(feature, label, w0, lam)
        while (m1 < 20):
            newf = get_result(feature, label, (w0 + rho ** m1 * dk), lam)
            oldf = get_result(feature, label, w0, lam)
            if newf < oldf + sigma * (rho ** m1) * (gk.T * dk)[0, 0]:
                mk = m1
                break
            m1 = m1 + 1
        w = w0 + rho ** mk * dk
        if k > m:
            s.pop(0)
            y.pop(0)
        sk = w - w0
        qk = get_gradient(feature, label, w, lam)
        yk = qk - gk
        s.append(sk)
        y.append(yk)

        t = len(s)
        a = []

        for i in range(t):
            alpha = (s[t - i -1].T * qk) / (y[t - i - 1].T * s[t - i - 1])
            qk = qk - alpha[0, 0] * y[t - i -1]
            a.append(alpha[0, 0])
        r = H0 * qk

        for i in range(t):
            beta = (y[i].T * r) / (y[i].T * s[i])
            r = r + s[i] * (a[t - i - 1] - beta[0, 0])
        if yk.T * sk > 0:
            dk = -r
        k = k + 1
        w0 = w
    return w0

           

中間還要有一個判斷下降方向的過程。

資料比較少，是以效果差别都不大。

總結

首先提到的是梯度下降，梯度下降算法雖然很簡單，但是下降的方向會有所偏差，可能胡局部不穩定，速度不會特别快，但是最終是會到達終點。于是改進一下，梯度下降是一階拟合，那麼換牛頓法二階拟合，但是牛頓法問題來了，疊代的方向有可能是錯誤的，是以改進一下，加點阻力，就算是不準确的，用linear search也可以調整一下。但是對于阻尼牛頓法還是存在了計算複雜度的問題，于是改進一下，用DFP直接做近似逆矩陣，達到了降低複雜度的問題，但是對于方向梯度的問題還是不是特别好，于是又出來了BFGS，用了比較牛逼的Sherman-Marrion公式求出了逆矩陣。問題又來了，每次都存儲這麼大的矩陣，有沒有辦法降低一些呢？于是改進了一下，就出現了LBFGS，用兩個nx1的矩陣來表示逆矩陣，并且隻存儲後m個更新的内容，改進一下出現了two-loop-recursion算法。後面又繼續改進了一下。這部分的知識比較靠近數值優化我心髒已經受不了了。

GitHub代碼： https://github.com/GreenArrow2017/MachineLearning/tree/master/MachineLearning/Linear%20Model/LogosticRegression/LinearRegaression