Optimization of Machine Learning

時間 2019-12-05

標籤 optimization machine learning 简体版

原文原文鏈接

機器學習就是須要找到模型的鞍點，也就是最優勢。由於模型不少時候並非徹底的凸函數，因此若是沒有好的優化方法可能會跑不到極值點，或者是局部極值，甚至是偏離。因此選擇一個良好的優化方法是相當重要的。首先是比較常規的優化方法：梯度降低。如下介紹的這些算法都不是用於當個算法，能夠試用於能可微的全部算法。 ###Gradient Descent 常見會用在logistics regression或者是linear regression裏面。好比logistics regression，這個模型直接等於0是求不出解析解的，全部只能用一些迭代方法求解。當樣本的label是 $\{1, 0\}$ 每個樣本正確機率： $P(y|x;\theta) = (h(x_i)^{y_i}(1-h(x_i))^{1-y_i})$ 最大似然函數化簡： $l(\theta) = \sum_{i = 0}^{m}(y_ilogh(x_i) + (1-y_i)logh(x_i))$ target function就是如上的 $l(\theta)$ ，要作的就是優化上面的損失函數。 $\frac{d(l(\theta))}{d\theta} = -\frac{1}{m}\sum_{i=0}^{m}\{y_i\frac{1}{h(x_i)}\frac{d}{d\theta}h(x_i)-(1-y_i)\frac{1}{1-h(x_i)}\frac{d}{d\theta}h(x_i)\}$ $= -\frac{1}{m}\sum_{i=0}^{m}\{y_i\frac{1}{g(\theta^Tx_i)}-(1-y_i)\frac{1}{1-g(\theta^Tx_i)}\}\frac{d}{d\theta}g(\theta^Tx_i)$ $= -\frac{1}{m}\sum_{i=0}^{m}\{y_i\frac{1}{g(\theta^Tx_i)}-(1-y_i)\frac{1}{1-g(\theta^Tx_i)}\}g(\theta^Tx_i)(1-g(\theta^Tx_i))\frac{d}{d\theta}\theta^Tx_i$ $=-\frac{1}{m}\sum_{i=0}^{m}\{y_i(1-g(\theta^Tx_i)) - (1-y_i)g(\theta^Tx_i)\}x_i^j$ $= \frac{1}{m}\sum_{i=0}^{m}(h(x_i)-y_i)x_i^j$ 這個就是當前點的梯度，而函數的更新方向就是向着梯度的負方向進行的，因此最後是要減：git

\theta_j = \theta_j - \alpha\frac{1}{m}\sum_{i=0}^{m}(h(x_i)-y_i)x_i^j

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

點，而在

點的梯度是上升的，因此負方向纔是正確的，最後獲得的梯度就是要減去了。這個就是gradient descent。對於步長 $\alpha$ ，能夠選擇 $\frac{\alpha}{ \sqrt{t}}$ ，t就是循環次數，隨着迭代次數的增多，愈來愈靠近極值點，再相同的迭代步長是有可能越過極值點的。事實上步長的選擇也是一個很重要的方法，後面會使用到。梯度降低是用給一個平面去擬合，也就是一次曲線去擬合，由於梯度降低是用泰勒一階展開的，因此是一次曲線。這樣形成的問題就是，gradient descent方法的視野很小，由於它只有一階擬合，因此當前找到那個方向並不必定是全局最優的方向，多是有點偏離，因此常常出現之字形的走勢，迭代的速度也不快。它只會看當前一階導數哪一個最適合，而不會去管二次或者三次的。 **梯度降低有三種：全量梯度降低，全部是數據一塊兒作迭代更新。數據量很大的話速度回很慢。小批量梯度降低，一部分一部分數據的就行迭代數據。隨機梯度降低SGD，隨機選取幾個進行迭代，可能迭代的方向會有誤差，可是隨着時間流逝大方向仍是同樣的。**代碼實現前面的logistics regression中已經有了，再也不重複。給予這種不足，給出了牛頓法的改進。 ###NewTon Method 既然梯度降低作的是一階擬合，那麼試一下二階擬合。泰勒展開二階： $f(x) = f(x_0) + f'(x_0)(x - x_0) + \frac{f''(x_0)}{2!}(x-x_0)^2$ 兩邊對

求導數：

$x - x_0 = -\frac{f'(x_0)}{f''(x_0)}$ $x = x_0 - \frac{f'(x_0)}{f''(x_0)}$ 代換一下： $x_{n+1} = x_n - \frac{f'(x_n)}{f''(x_n)}$ 牛頓法有點像座標上升的方法，先迭代一個座標，再迭代另外一個，和SMO相似。能夠看到牛頓法這個時候就是二次曲線的擬合了：

左邊的就是成功擬合的，右邊的那個就是方向錯誤，擬合沒有成功的例子。因此牛頓法不必定會按照正確的方向擬合。上面牛頓法的式子是對於單變量的，若是是對變量，那麼下面的二階導要用到Hession矩陣。因此對於多變量的牛頓法：

$d_k = -H_k^{-1}g_k$ 稱爲是牛頓方向，剛剛提到，迭代方向要和一階導數方向相反，也就是

，這個想法和剛剛提出的沿着梯度方向的反方向降低的思路是同樣的。分解一下這個公式： $g_k^TH_k^{-1}g_k > 0$ ，因此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函數，使用牛頓法代替梯度降低作邏輯迴歸。github

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
複製代碼

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

'''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])
複製代碼

運行函數。app

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)

複製代碼

生成數據。dom

若是是隨機值的話，多試幾回有可能方向是錯誤的。這裏能跑這麼快是由於初始值設置的好。對於牛頓法的這個問題仍是有改進方法的。對於 $\alpha$ 步長的算法研究是能夠解決這個問題。既然有時候牛頓法降低的方向或許是不對的，可使用一個阻力來阻擋一下，減慢速度，等到下一次能夠從新選擇。因此接下來就是對於 $\alpha$ 步長的研究。 ###線搜索步長最小化意味着 $f(\alpha) = minimizef(x_k + \alpha d_k)$ ，

就是方向了,這裏的

都是常數，因此能夠看作是一個 $\alpha$ 的函數。 #####二分搜索算法

$\alpha = argmin_{\lambda>0}f(x+\lambda p)$ 假設 $g(\alpha)=f(x_k + \alpha d_k)$ ,若是是一個凸函數，須要 $g'(\alpha) = 0$ ，由於 $g'(\lambda) = \nabla f(x_k + \lambda p)p$ ，p就當前的梯度降低的方向，因此向着梯度的反方向得：

，假如咱們知道一個 $g'(\lambda)>0$ ，那麼就能夠知道確定有一個

使得 $g'(x_k) = 0,x_k \in (0, x_k)$ ，因此二分搜索的流程很簡單： ①k = 0， $\alpha_0 = 0$ ， $\alpha_1 = a$ ② $m = \frac{\alpha_0 + \alpha_1}{2}$ ，若是 $g'(m) > 0,\alpha_1 = m,k+=1$ 若是 $g'(m) < 0,\alpha_0 = m,k+=1$ 若是

，不然就繼續循環。二分搜索也叫精確搜索，由於這樣不斷的迭代下去是能夠很精確的搜索到一個合適的值，使用二分查找法來求步長的計算複雜度很高, 由於在最小化f(x)f(x)的每次迭代中咱們都須要執行一次線搜索, 而每次線搜索都要用上述的二分查找算法，由於這個精確的值在一些凸函數裏面會是很小很小的，甚至後面多加一個小數，這樣算起來複雜度很是高的。咱們能夠在犧牲必定的精度的條件下來加快計算速度, 回溯線搜索是一種近似線搜索算法。 #####回溯線搜索對於上訴二分法的問題，這裏提出了改進，犧牲掉一些精度，求一個大概的範圍就好，因而有：

f(x_k+\alpha p_k) <= f(x_k) + c_1 \alpha \nabla f_k^Tp_k$$上訴條件稱爲**充分降低條件**，$c_1 \in (0, 1)$
$$f(x_k+\alpha p_k) >= f(x_k) + c_1 (1-\alpha) \nabla f_k^Tp_k$$這兩個準則就是Armijo-Goldstein準則，①目標函數值應該有足夠的降低；②一維搜索的步長α不該該過小。這兩個目標其實很明顯，足夠的降低其實就是式子1：![](https://upload-images.jianshu.io/upload_images/10624272-4e8a494e2452c6e8.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)能夠看到接受這個條件有兩個區間，有時候會選擇到第一個區間的內容，也就是第一個區間的內容，因此第二條式子的做用就是捨得步長不要過小了。
>>####爲何$c_1 \in (0, \frac{1}{2})$
>> **①對於$c_1$的選擇是必定要在(0,0.5)之間的，若是沒有這個條件，會影響算法的超線性收斂性。這個怎麼收斂的我就不知道了，還得看看其餘的凸優化書籍。
②$f(x_k + \alpha p_k)$作泰勒展開，$f(x_k + \alpha p_k) = f(x_k) + \alpha g_k^T d_k+o(\alpha)$，去掉高階無窮小，$f(x_k + \alpha p_k) = f(x_k) + \alpha g_k^T d_k$，和上式就是相差了一個$p$而已，這樣就保證了小於的這個條件。
③因爲$p \in (0, \frac{1}{2})$，因此$(1-p) \in (\frac{1}{2}, 1)$，並且$g_k^Td_k<0$，因此有$\alpha_k(1-p)g_k^Td_k < \alpha_k p g_k^Td_k < 0$，這就保證了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

```
求海塞矩陣換了一個寫法，好看一點，可是效果都是同樣的。
線性迴歸的損失函數：$$loss = \frac{1}{2}\sum_{i=1}^{m}(y^i-\sum_{j-0}^{n-1}w_j*x_j^i)^2

線性迴歸的一階導數：$$\frac{dloss}{dw_j} = -\sum_{i=1}^{m}(y^i - \sum_{j=0}^{n-1}w_j*x_j^i)x_j^i$$ 線性迴歸的二階導數：$$\frac{dloss}{dw_jdw_k} = \sum_{i=1}^m(x_j^ix_k^i)$$ 對照上面的公式是能夠一一對應的。對於Armrji搜索只是用了第一條公式：$$f(x_k+\alpha p_k) <= f(x_k) + c_1 \alpha \nabla f_k^Tp_k$$一旦超過這個限制就能夠退出了。機器學習

效果其實仍是能夠的。固然，也能夠三條準則都使用上，只不過使用一條通常也足夠了。到這裏線搜索尚未完。 #####Wolfe-Powell準則有Armijo-Goldstein準則仍是不夠的，以下圖：

能夠看到(a,b),(c,d)就是選擇的區間，可是很明顯這些區間已經避開了最低的點，固然這不是必定會，可是有這個可能，爲了解決這個問題就出現了Wolfe-Powell準則。 Wolfe-Powell準則也有兩個表達式，第一個表達式就是$$f(x_k+\alpha p_k) <= f(x_k) + c_1 \alpha \nabla f_k^Tp_k$$第二條式子就去掉了，由於它會很容易的把極值點排除在外，因此並不須要這個條件。添加的就是：$$\nabla f(x_k + \alpha_k d_k)^Td_k >= c_2 g_k^Td_k,c_2 \in (c_1, 1)$$這條式子的直觀解釋就是， 當前可選點的斜率要大於初始點斜率的倍。，而初始點（ α=0 的點）處的切線是比 e 點處的切線要「斜」的，因爲

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

e和f後面的點就是知足的了， $f(x_k+\alpha p_k) >= f(x_k) + c_1 (1-\alpha) \nabla f_k^Tp_k$ 這個條件是沒有的了。這個條件仍是不太好，有時候會給出更嚴謹的一個條件： $|\nabla f(x_k + \alpha_k d_k)^Td_k| >= -c_2 g_k^Td_k,c_2 \in (c_1, 1)$ ，其實就是兩邊都給絕對值，可是右邊的原本就是一個負數了，加上符號就是正數。這個條件就更強了：

直接限定在了(e,g)和(f,h)範圍以內了。這就是整個Wolfe-Powell準則。 #####線搜索的總結線搜索到這裏基本就結束了。尋找一個最佳步長就是要求優化以後的函數值是要最小的，因而給出的目標函數就是： $\alpha = argmin_{\lambda>0}f(x+\lambda p)$ 。解決這個問題有兩種方法，一種就是二分法的查找，這種方法可能很是的精確搜索到一個最佳的值，可是他的計算複雜度有點高；有時候咱們其實不太須要太過精確的值，咱們只是須要一個大概模糊的就行了，因而出現了回溯搜索，這是屬於模糊搜索，給出的就是Armijo-Goldstein準則，可是這兩個準則也有很差的地方，有時候會把極值點排除在外，因而引入曲率條件，把Armijo-Goldstein準則的第二條式子換掉，就出現了Wolfe-Powell準則，條件再強一點兩邊加上絕對值，這樣就出現了最終的Wolfe-Powell準則。對於步長的研究並不針對某一個算法，對於優化降低的算法均可以使用，梯度降低，牛頓法，擬牛頓法均可以用到，具備很強的普適性。梯度降低的步長也是能夠經過這種方式進行選擇最優步長，牛頓法用Armijo搜索的方法是能夠獲得全局牛頓法，也叫阻尼牛頓法，這樣可使得迭代方向能夠避免向錯誤的方向進行，增長點阻力。這篇文章主要的研究對象仍是牛頓法，因此下面的三個算法分別就是DFP，BFGS，LBFGS。 ###DFP 前面的阻尼牛頓法解決的就是對於迭代方向的問題，可是對於計算複雜度這個問題尚未獲得解決，由於矩陣求擬是一個很大的工程量，若是維度一可能是計算複雜度是很大的，因此擬牛頓法基本上都是構造一個和Hession矩陣類似的矩陣來替代。可是問題來了，用近似矩陣替代，效果可能會很差，給出以下證實：首先由牛頓法能夠獲得： $-\frac{\nabla f(x_i)}{H} = x - x_i$ $\nabla f(x_i) = -(x-x_i)H$ $\nabla f(x_i)d = (x-x_i)H * H^{-1}g$ $\nabla f(x_i)d = -(x-x_i)H(x-x_i)$ 因爲 $\nabla f(x_i) d < 0$ ，因此右邊也是要求小於0，那麼天然就是大於0了，這也就證實了Hession矩陣的正定性。在逐步尋找最優解的過程當中要求函數值是降低的，可是有時候Hession矩陣不必定是半正定的，可能會使得函數值不降反升。擬牛頓法可使目標函數值沿降低方向走下去，而且到了最後，在極小值點附近，可以使構造出來的矩陣與Hesse矩陣「很像」了，這樣，擬牛頓法也會具備牛頓法的二階收斂性。 #####推導證實涉及到Hession矩陣，天然就涉及到二階導數矩陣了，Hession矩陣只是對於多變量二階導數的一種描述方式。Taylor展開：$$f(x) = f(x_{i+1})+(x-x_{i+1})^T\nabla f(x_{i+1}) + \frac{1}{2}(x-x_{i+1})^TH_{i+1}(x+ x_{i+1})+o(x+x_{i+1})....$$仍是用二次導數作近似，不少函數均可以用二次導數作近似。去掉高次項兩邊求導數：$$\nabla f(x_i) = \nabla f(x_{i+1}) + H_{i+1}(x_i-x_{i+1}) \ H^{-1}\nabla f(x_i)=H^{-1}\nabla f(x_{i+1}) + (x_i-x_{i+1}) \ H^{-1}(\nabla f(x_i) - \nabla f(x_{i+1})) = (x_i - x_{i+1})$$ 上面那個方程就是擬牛頓方程。上面那個矩陣就是一個近似矩陣。近似矩陣有不少種，若是按照上面那種方式進行計算，複雜度沒有降多少，因此比較常見的方法就是迭代計算，比較常見的迭代方式：$$H_{i+1} = H_i + E_i$$ 那麼接下來問題來了，

要這麼設計。通常給出的設計就是$$E = mvv^T+nww^T$$由於但願每次迭代h能有一個微小的變化而不是巨大的變化，這樣纔有可能收斂。並且這個結構設計的也很簡單，也是對稱的。代入上面的擬牛頓方程：$$H_{i+1}(\nabla f(x_{i+1})-\nabla f(x_i)) = (H_i+E_i)(\nabla f(x_{i+1}) - \nabla f(x_i)) \ (H_i + mvv^T+nww^T)(\nabla f(x_{i+1}) - \nabla f(x_i))=x_{i+1} - x_i \ H_i(\nabla f(x_{i+1})-\nabla f(x_i)) + mvv^T(\nabla f(x_{i+1})-\nabla f(x_i)) + nww^T(\nabla f(x_{i+1})-\nabla f(x_i)) = x_{i+1}-x_i \ support: q_i = (\nabla f(x_{i+1})-\nabla f(x_i)) ,s_i = x_{i+1}-x_i \ then: H_iq_i+mvv^Tq_i+nww^Tq_i = s_i \ H_iq_i+v(mv^Tq_i)+w(nw^T)q_i = s_i$$仔細看一下這個式子，右邊的

是一個nx1的向量，而

均爲實數，也就是點積。能夠直接假設

，因而有

。代入上面的式子：$$H_{i+1} = H_i+\frac{s_is_i^T}{s_i^Tq_i}-\frac{H_iq_iq_i^TH_i}{q_i^TH_iq_i}$$一開始我看到這個推導過程，我是一臉懵逼的。原來數學家也有猜的時候，這個過程和前面假設H的形式在凸優化裏面，是用"很顯然"來描述的，嗯，很顯然！我就是他媽的沒看出來顯然在哪了，後面的LBFGS也是，真的很顯然，可是這個我是真的不知道。看到網上的不少解釋都離不開特殊狀況，也就是爲0或者是直接假設的，都差很少的解釋。已知初始正定矩陣H，從一個初始點開始（迭代），用式子

來計算出下一個搜索方向，並在該方向上求出可以使目標函數極小化的步長α，而後用這個步長，將當前點挪到下一個點上，並檢測是否達到了程序停止的條件，若是沒有達到，則用上面所說的[13]式的方法計算出下一個修正矩陣H，並計算下一個搜索方向……周而復始，直到達到程序停止條件。因此整一個流程： ①給定一個初值 ②搜索方向： ③利用搜索獲得步長 $\alpha$ ，可使用上面提到的Armrji搜索或者等等的改進方法。 ④更新一波 ⑤計算 $y_k = g_{k+1}-g_{k}$ ⑥ $H_{i+1} = H_i+\frac{s_is_i^T}{s_i^Tq_i}-\frac{H_iq_iq_i^TH_i}{q_i^TH_iq_i}$ ，轉回去繼續更新。 DFP仍是有缺點的：

具體主要內容就是說若是Hession矩陣的是錯誤估計的，大家BFGS會在不多的迭代回到正確的方向，可是DFP在這方面的效果不明顯。至於爲何不明顯，我從公式還不能直接看出來。既然提到了BFGS，下面就是BFGS了。 ###BFGS BFGS和DFP實際上是屬於對偶的解法，一個直接求海賽矩陣逆矩陣（DFP），另外一個就是求海賽矩陣（BFGS）。仍是使用擬牛頓公式：

，而DFP使用的是 $H_i^{-1}q_i = s_i$ ，直接求出來的就是逆矩陣了。推導過程其實很簡單，和上面DFP相似，差異不大，甚至除了換一下沒什麼差異的。 #####推導證實假設都和上面DFP的同樣。仍是 $H_{i+1} = H_{i} + E_i$ ，E就是矯正矩陣。E的形式和以前的是同樣的

。代入上面的公式：$$q_i = H_{i+1}s_i \ q_i = (H_i + mvv^T + nww^T)s_i \ q_i = H_is_i+(mv^Ts_i)v+(nw^Ts_i)w$$一樣假設

，一樣選取特殊狀況， $v = q_i,w = H_is_i,m = \frac{1}{v^Ts_i},n = -\frac{1}{w^Ts_i}$ ，代進矯正函數的表達式：$$H_{i+1} = H_i+\frac{vv^T}{v^Ts_i} - \frac{ww^T}{w^Ts_i} \ H_{i+1} = H_i+\frac{q_iq_i^T}{q_i^Ts_i} - \frac{H_is_is_i^TH_i^T}{s_i^TH_is_i} $$仔細看一些這個式子，和DFP那個式子是換了一下位置而已，因此這個過程和步驟真的沒有什麼好說的。按照常規套路，通常是這樣：因此整一個流程： ①給定一個初值 ②搜索方向： $d_k = -H_k^{-1} *g_k$ ③利用搜索獲得步長 $\alpha$ ，可使用上面提到的Armrji搜索或者等等的改進方法。 ④更新一波 ⑤計算 $y_k = g_{k+1}-g_{k}$ ⑥ $H_{i+1} = H_i+\frac{q_iq_i^T}{q_i^Ts_i} - \frac{H_is_is_i^TH_i^T}{s_i^TH_is_i}$ ，轉回去繼續更新。然而，若是是這樣，複雜度仍是存在的，仍是得求個導數啊。因此這樣的作法天然是不行的，由於擬牛頓法的改進有一個很重要的誘因就是計算複雜度的問題了。 對於求逆矩陣，有一個很是重要的公式——Sherman-Morrison公式：$$(A+uv^T)^{-1}=A^{-1}-\frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}$$ 用Sherman-Morrison公式改造一下： $H_{k+1}^{-1}=(I-\frac{s_kq_k^T}{q_k^Ts_k})H_k^{-1}(I-\frac{q_ks_k^T}{q_k^Ts_k})+\frac{s_ks_k^T}{q_k^Ts_k}$ 至因而怎麼改造的，有興趣本身看吧，原諒我並無看懂，懂了我也想不出來。根據這個結果改造一下：因此整一個流程： ①給定一個初值 ②搜索方向： $d_k = -H_k^{-1} *g_k$ ③利用搜索獲得步長 $\alpha$ ，可使用上面提到的Armrji搜索或者等等的改進方法。 ④更新一波 ⑤計算 $y_k = g_{k+1}-g_{k}$ ⑥ $H_{k+1}^{-1}=(I-\frac{s_kq_k^T}{q_k^Ts_k})H_k^{-1}(I-\frac{q_ks_k^T}{q_k^Ts_k})+\frac{s_ks_k^T}{q_k^Ts_k}$ ，轉回去繼續更新。完美的達到了下降複雜度的目標。可是我仍是不知道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算法作了近似，而是存儲計算時所須要的 $\{s_k\},\{q_k\}$ ，而後利用哪兩個向量的計算來近似代替。對於 $\{s_k\},\{q_k\}$ 也不是存儲全部的，而是隻保存那麼後面的m個，也就是最新的m個就好。這樣就把空間複雜度降到。毫無疑問，LBFGS是從BFGS演化而來的。對於逆矩陣的公式：$$H_{k+1}^{-1}=(I-\frac{s_ky_k^T}{y_k^Ts_k})H_k^{-1}(I-\frac{y_ks_k^T}{y_k^Ts_k})+\frac{s_ks_k^T}{y_k^Ts_k}$$ 假設 $p_k = \frac{1}{y_k^Ts_k},v_k = 1-p_ky_ks_k^T$ ，改造一下上式：$$H_{k+1}^{-1} = v_k^TH_k^{-1}v_k+p_ks_ks_k^T$$ 假設： $D_k = H_k^{-1}$ 遞推一下：oop

D_1 = v_0^TD_0v_0+p_0s_0s_0^T\\
D_2 = v_1^TD_1v_1+p_1s_1s_1^T\\
= v_1^T(v_0^TD_0v_0+p_0s_0s_0^T)v_1+p_1s_1s_1^T\\
= v_1^Tv_0^TD_0v_0v_1+v_1s_0s_0^Tv_1+p_1s_1s_1^T

D_3 = v_2^Tv_1^Tv_0^TD_0v_0v_1v_2+v_2^Tv_1^Tp_0s_0s_0^Tv_1v_2+v_2^Tp_1s_1s_1^Tv_2+p_2s_2s_2^T

根據推論，通常的有：$$D_{k+1} = (v_k^Tv_{k-1}^T...v_1^Tv_0^T)D_0(v_0v_1...v_{k-1})\ +(v_k^Tv_{k-1}^T...v_2^Tv_1^T)(p_0s_0s_0^T)(v_1v_2...v_{k-1}v_k)\ +(v_k^Tv_{k-1}^T...v_3^Tv_2^T)(p_1s_1s_1^T)(v_2v_3...v_k-1v_k)\ +...\ +(v_k^Tv_{k-1}^T)(p_{k-2}s_{k-2}s_{k-2}^T)(v_{k-1}v_{k})\ +v_k^T(p_{k-1}s_{k-1}s_{k-1}^T)v_k\ +p_ks_ks_k^T$$ 能夠看的出，當前的海賽矩陣的逆矩陣是能夠由 $\{s_i,q_i\}$ 給出的，因此直接存儲對應的便可，上面提到過爲了壓縮內存，能夠只是存儲後面的m個，因此更新一下公式：$$D_{k+1} = (v_k^Tv_{k-1}^T...v_{k-m+1}^Tv_{k-m}^T)D_0(v_{k-m}v_{k-m+1}...v_{k-1})\ +(v_k^Tv_{k-1}^T...v_{k-m+2}^Tv_{k-m+1}^T)(p_0s_0s_0^T)(v_{k-m+1}v_{k-m+2}...v_{k-1}v_k)\ +(v_k^Tv_{k-1}^T...v_{k-m+3}^Tv_{k-m+2}^T)(p_1s_1s_1^T)(v_{k-m+2}v_{k-m+3}...v_k-1v_k)\ +...\ +(v_k^Tv_{k-1}^T)(p_{k-2}s_{k-2}s_{k-2}^T)(v_{k-1}v_{k})\ +v_k^T(p_{k-1}s_{k-1}s_{k-1}^T)v_k\ +p_ks_ks_k^T$$ 公式長是有點長，可是已經有大神給出了一個很好的算法解決這個問題，two-loop-recursion算法： $y_k=\nabla f(x_k)-\nabla f(x_{k-1})$ $q = \nabla f(x_k)$ $\alpha_i = p_{k-i}s_{k-i}^Tq$ $q = q-\alpha_iy_{k-i}$ end for $r = H_{k}^0q$ $\beta = p_{k-1} y_{k-1}^T r$ $r = r+s_{k-i}(\alpha_i - \beta)$ 這個式子到底行不行呢？證實一下理論：學習

這樣就證實這個算法的正確性。然而其實我根本不關心這個算法正確性，我只是想知道 ###這是怎麼想出來的，說實話第一眼看根本沒有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代碼：github.com/GreenArrow2…

相關標籤/搜索

每日一句

每一个你不满意的现在，都有一个你没有努力的曾经。