深度學習中的優化器比較

一. 幾個數學概念

1) 梯度（一階導數）

考慮一座在 (x1, x2) 點高度是 f(x1, x2) 的山。那么，某一點的梯度方向是在該點坡度最陡的方向，而梯度的大小告訴我們坡度到底有多陡。

2) Hesse 矩陣（二階導數）

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

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

其中 x為列向量，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 矩陣（一階偏導與二階偏導的關系）。

二.優化算法介紹

補充在前：實際上在我使用LSTM為流量基線建模時候，發現有效的激活函數是elu、relu、linear、prelu、leaky_relu、softplus，對應的梯度算法是adam、mom、rmsprop、sgd，效果最好的組合是：prelu+rmsprop。

1.批量梯度下降（Batch gradient descent，BGD）

θ=θ−η⋅∇θJ(θ) 

每迭代一步，都要用到訓練集的所有數據，每次計算出來的梯度求平均 

η代表學習率LR

2.隨機梯度下降（Stochastic Gradient Descent，SGD）

θ=θ−η⋅∇θJ(θ;x(i);y(i)) 

通過每個樣本來迭代更新一次，以損失很小的一部分精確度和增加一定數量的迭代次數為代價，換取了總體的優化效率的提升。增加的迭代次數遠遠小於樣本的數量。

缺點：

• 對於參數比較敏感，需要注意參數的初始化 
• 容易陷入局部極小值 
• 當數據較多時，訓練時間長 
• 每迭代一步，都要用到訓練集所有的數據。

import tensorflow as tf
tf.compat.v1.train.GradientDescentOptimizer()

3. 小批量梯度下降（Mini Batch Gradient Descent，MBGD）

θ=θ−η⋅∇θJ(θ;x(i:i+n);y(i:i+n)) 

為了避免SGD和標准梯度下降中存在的問題，對每個批次中的n個訓練樣本，這種方法只執行一次更新。【每次更新全部梯度的平均值】

在深度學習中為了解決學習率的靈活性，防止迭代到一定次數時，因學習率過大使得函數值在某個值的附近來回波動而不能繼續先最優解的方向迭代，可以對學習率進行指數退化。

import tensorflow as tf
'''
tf.train.exponential_decay
#該函數返回以下結果
decayed_learning_rate = learning_rate *
         decay_rate ^ (global_step / decay_steps)
##例： 以0.96為基數，每100000 步進行一次學習率的衰退
'''
global_step = tf.Variable(0, trainable=False)
starter_learning_rate = 0.1
learning_rate = tf.train.exponential_decay(starter_learning_rate, global_step,
                                           100000, 0.96, staircase=True)
# Passing global_step to minimize() will increment it at each step.
learning_step = (
    tf.train.GradientDescentOptimizer(learning_rate)
    .minimize(...my loss..., global_step=global_step)
)

4.近端梯度下降法 (proximal gradient descent)

近端梯度下降方法，這個是tf取的名字，其實在文章中，作者簡稱該梯度下降算法為folos，之所以叫接近的梯度下降算法，是因為作者在原梯度下降算法的基礎上加入了修正，進而解決了一些不可微的問題，

近端梯度下降法是眾多梯度下降 (gradient descent) 方法中的一種，其英文名稱為proximal gradident descent，其中，術語中的proximal一詞比較耐人尋味，將proximal翻譯成“近端”主要想表達"（物理上的）接近"。與經典的梯度下降法和隨機梯度下降法相比，近端梯度下降法的適用范圍相對狹窄。對於凸優化問題，當其目標函數存在不可微部分（例如目標函數中有  -范數或跡范數）時，近端梯度下降法才會派上用場。

特定的凸優化問題

一般而言，近端梯度下降法常用於解決以下這類優化問題：

假設目標函數  是由  和  疊加而成，其中，限定  是可微的凸函數、  是不可微 (或局部不可微) 的凸函數。

以線性回歸問題為例

給定  ，若線性回歸表達式為  ，其中，變量  表示系數向量。優化的目標函數可以寫成如下形式：

其中，  表示  -范數；  表示  -范數；正則項中的  -范數是為了保證變量 是稀疏的。在這個目標函數中，正則項中的  -范數便是一個不可微的凸函數。由此便引申出采用近端梯度下降法對這類優化問題進行求解。

近端算子 (proximal operator)

一般來說，對於任意向量  ，如果不可微函數為，則近端算子 (proximal operator) 為

其中，  表示關於變量  和函數  的近端算子； 表示關於變量 軟閾值函數（下面將着重介紹），選用符號  是因為軟閾值soft-thresholding的英文首字母為s。

這條公式想表達的意思是：對於任意給定的  ，我們希望找到使得  最小化的解  ，其中，  是一個新增參數，表示近端梯度下降的步長 (step size)。實際上，近端算子只和不可微的凸函數  有關 (如公式(2)中的)。

當不可微的凸函數的形式為  時，則對應的軟閾值函數為

如公式(2)，當  （  -范數作為正則項，其正則項參數為  ）時，軟閾值函數則為，其定義為

公式(4)給出的軟閾值函數恰好與公式(2)中的目標函數相對應。至於為什么軟閾值函數的形式是這樣，其推導則需自己百度查找。

近端梯度下降 (proximal gradient descent)

對於優化問題  ，變量  的迭代遞推公式為

其中，變量上標的  表示當前迭代次數。

迭代遞推公式證明過程（涉及知識：泰勒展開）:

注意：由於公式第三行中的  和第四行中的  均與決策變量  無關，因此公式第三行等於公式第四行。

線性回歸問題的迭代過程

回到公式(1)，我們定義了函數  的表達式為

由於  ，根據近端梯度下降法的迭代遞推公式可以推導出變量的近端梯度更新公式為

這條更新公式通常也被稱為迭代軟閾值算法 (iterative soft-thresholding algorithm, ISTA)。

眾所周知，如果將這里的  換成  ，則整個線性回歸的優化問題就會變成一個最小二乘問題。獲取變量  的最優解也會變得非常簡單，只需要令  即可，此時的最優解為

其中，  表示正則項參數。

import tensorflow as tf
tf.compat.v1.train.ProximalGradientDescentOptimizer()

當然，這里還有一個很值得思考的問題：

• 對於上述線性回歸問題，為何選用近端梯度下降法，而非次梯度 (subgradient descent) 下降法呢？

實際上，次梯度下降法也能求解這種線性回歸問題。一般而言，對於損失函數

我們可以得到正則項的梯度為  然而不幸的是，當  時，正則項不可微。為了消除這個影響，我們可以采用次梯度  ，當再遇到  時，次梯度就會等於  .

相應地，損失函數的次梯度為

與梯度下降法類似，我們可以用次梯度下降法使得損失函數最小化。但仔細看一下損失函數次梯度的表達式  ，我們會發現：次梯度法可能不會像我們所期許的那樣，得到真正的稀疏解 ，次梯度法只會使解的部分數值 (即  的部分元素) 盡可能小 (盡可能接近  )。因此，只有選用近端梯度下降這類方法，我們才能獲得稀疏解，並同時提升迭代過程的收斂性。

5.指數加權平均的概念

指數加權平均(exponentially weighted averges)，也叫指數加權移動平均，是一種常用的序列數據處理方式。

它的計算公式如下：

其中，

• θ_t：為第 t 天的實際觀察值，
• V_t: 是要代替 θ_t 的估計值，也就是第 t 天的指數加權平均值，
• β： 為 V_{t-1} 的權重，是可調節的超參。( 0 < β < 1 )

直接看上面的數據圖會發現噪音很多，

這時，我們可以用 指數加權平均 來提取這組數據的趨勢。

 

可以看出，紅色的數據比藍色的原數據更加平滑，少了很多噪音，並且刻畫了原數據的趨勢。

指數加權平均，作為原數據的估計值，不僅可以 1. 撫平短期波動，起到了平滑的作用，2. 還能夠將長線趨勢或周期趨勢顯現出來。

所以應用比較廣泛，在處理統計數據時，在股價等時間序列數據中，CTR 預估中，美團外賣的收入監控報警系統中的 hot-winter 異常點平滑，深度學習的優化算法中都有應用。

import matplotlib.pyplot as plt
import numpy as np
# from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits import mplot3d
np.random.seed(28)
#生成每個月的平均氣溫數據
average_temperature = [13,10,13,18,23,28,33,36,33,28,]
temperature1 = np.random.normal(loc=13,scale=3,size=(31,))
temperature2 = np.random.normal(loc=10,scale=3,size=(28,))
temperature3 = np.random.normal(loc=13,scale=3,size=(31,))
temperature4 = np.random.normal(loc=18,scale=3,size=(30,))
temperature5 = np.random.normal(loc=23,scale=3,size=(31,))
temperature6 = np.random.normal(loc=28,scale=3,size=(30,))
temperature7 = np.random.normal(loc=33,scale=3,size=(31,))
temperature8 = np.random.normal(loc=36,scale=3,size=(31,))
temperature9 = np.random.normal(loc=33,scale=3,size=(30,))
temperature10 = np.random.normal(loc=28,scale=3,size=(31,))
temperature11 = np.random.normal(loc=23,scale=3,size=(30,))
temperature12 = np.random.normal(loc=18,scale=3,size=(31,))
temperature_data = np.hstack((temperature1,temperature2,temperature3,temperature4,temperature5,temperature6,temperature7,
                  temperature8,temperature9,temperature10,temperature11,temperature12))
days = np.arange(1,len(temperature_data)+1)

#求指數加權平均數
def exp_mov_avg(data,beta,amend=False):
    '''
    :param data: 一維數據
    :param beta: 加權系數
    :param amend: 是否修正
    :return:
    '''
    v=[]
    v_t = 0
    if amend:
        for i,e in enumerate(data):
            print(i)
            v_t = beta * v_t + (1 - beta) * e
            #修正數值
            amend_v = v_t/(1-beta**(i+1))
            v.append(amend_v)
    else:
        for e in data:
            v_t = beta*v_t+(1-beta)*e
            v.append(v_t)
    return v

mov_avg = exp_mov_avg(temperature_data,0.5)
mov_avg2 = exp_mov_avg(temperature_data,0.9)
mov_avg3 = exp_mov_avg(temperature_data,0.98)
mov_avg4 = exp_mov_avg(temperature_data,0.9,True)
plt.scatter(days,temperature_data)
plt.plot(days,mov_avg,'y-',label='beta=0.5')
plt.plot(days,mov_avg2,'r-',label='beta=0.9')
plt.plot(days,mov_avg3,'g-',label='beta=0.98')
plt.plot(days,mov_avg4,color='purple',ls='-',label='beta=0.9 and amend')
plt.legend()
plt.show()

'''使用tensorflow的api'''
# import tensorflow as tf
# tf.compat.v1.enable_eager_execution()
# var = tf.Variable(0,dtype=tf.float32)
# ema=tf.compat.v1.train.ExponentialMovingAverage(0.9)
# mov_avg = []
# for e in temperature_data:
#     new_var = tf.constant(e,dtype=tf.float32)
#     var.assign(new_var)
#     ema.apply([var])
#     mov_avg.append(ema.average(var).numpy())
# 
# plt.scatter(days,temperature_data)
# plt.plot(days,mov_avg,'r-',label='beta=0.9')
# plt.legend()
# plt.show()

β 如何選擇

根據前面的計算式子：

將 V_{100} 展開得到：

 
從這里我們就可已看出指數加權平均的名稱由來，第100個數據其實是前99個數據加權和，而前面每一個數的權重呈現指數衰減，即越靠前的數據對當前結果的影響較小 

這里可以看出，V_t 是對每天溫度的加權平均，之所以稱之為指數加權，是因為加權系數是隨着時間以指數形式遞減的，時間越靠近，權重越大，越靠前，權重越小。

再來看下面三種情況：

當 β = 0.9 時，指數加權平均最后的結果如圖紅色線所示，代表的是最近 10 天的平均溫度值；

當 β = 0.98 時，指結果如圖綠色線所示，代表的是最近 50 天的平均溫度值；

當 β = 0.5 時，結果如下圖黃色線所示，代表的是最近 2 天的平均溫度值；

β 越小，噪音越多，雖然能夠很快的適應溫度的變化，但是更容易出現奇異值。

β 越大，得到的曲線越平坦，因為多平均了幾天的溫度，這個曲線的波動更小。但有個缺點是，因為只有 0.02 的權重給了當天的值，而之前的數值權重占了 0.98 ，曲線進一步右移，在溫度變化時就會適應地更緩慢一些，會出現一定延遲。

通過上面的內容可知，β 也是一個很重要的超參數，不同的值有不同的效果，需要調節來達到最佳效果，一般 0.9 的效果就很好。

缺點：存在開始數據的過低問題，可以通過偏差修正，但是在深度學習的優化算法中一般會忽略這個問題 
 
當t不斷增大時，分母逐漸接近1，影響就會逐漸減小了

優點：【相較於滑動窗口平均】 

• 1.占用內存小，每次覆蓋即可 
• 2.運算簡單

6.Momentum（動量梯度下降法）

上面提到了一些 指數加權平均 的應用，這里我們着重看一下在優化算法中的作用。

動量優化算法在梯度下降法的基礎上進行改變，具有加速梯度下降的作用。類別：（1）標准動量優化方法Momentum,(2)NAG動量優化方法（NAG在Tensorflow中與Momentum合並在同一函數tf.train.MomentumOptimizer中，可以通過參數配置啟用。）

以 Momentum 梯度下降法為例，

Momentum 梯度下降法，就是計算了梯度的指數加權平均數，並以此來更新權重，它的運行速度幾乎總是快於標准的梯度下降算法。

這是為什么呢？

讓我們來看一下這個圖，

例如這就是我們要優化的成本函數的形狀，圖中紅點就代表我們要達到的最小值的位置，假設我們從左下角這里出發開始用梯度下降法，那么藍色曲線就是一步一步迭代，一步一步向最小值靠近的軌跡。可以看出這種上下波動，減慢了梯度下降法的速度，而且無法使用更大的學習率，因為如果用較大的學習率，可能會偏離函數的范圍。如果有一種方法，可以使得在縱軸上，學習得慢一點，減少這些擺動，但是在橫軸上，學習得快一些，快速地從左向右移移向紅點最小值，那么訓練的速度就可以加快很多。

這個方法就是動量 Momentum 梯度下降法，它在每次計算梯度的迭代中，對 dw 和 db 使用了指數加權平均法的思想。

momentum是模擬物理里動量的概念，積累之前的動量來替代真正的梯度。引入一個累計歷史梯度信息動量加速SGD。公式如下： 
 

然而網上更多的是另外一種版本，即去掉（1-β） 

 

相當於上一版本上本次梯度的影響權值*1/(1-β) 

兩者效果相當，只不過會影響一些最優學習率的選取 

優點

• 下降初期時，使用上一次參數更新，下降方向一致，乘上較大的μ能夠進行很好的加速
• 下降中后期時，在局部最小值來回震盪的時候，gradient→0，β得更新幅度增大，跳出陷阱
• 在梯度改變方向的時候，μ能夠減少更新

即在正確梯度方向上加速，並且抑制波動方向張的波動大小，在后期本次計算出來的梯度會很小，以至於無法跳出局部極值，Momentum方法也可以幫助跳出局部極值參數設置 。加上動量項就像從山頂滾下一個球，求往下滾的時候累積了前面的動量(動量不斷增加)，因此速度變得越來越快，直到到達終點。同理，在更新模型參數時，對於那些當前的梯度方向與上一次梯度方向相同的參數，那么進行加強，即這些方向上更快了；對於那些當前的梯度方向與上一次梯度方向不同的參數，那么進行削減，即這些方向上減慢了。因此可以獲得更快的收斂速度與減少振盪。

β的常用值為0.9，即可以一定意義上理解為平均了前10/9次的梯度。 至於LR學習率的設置，后面所有方法一起總結吧。

from matplotlib import pyplot as plt
import numpy as np

fig = plt.figure()
x = np.arange(-0.8, 1.2, 0.025)
plt.plot(x,-x**3-x**2+2*x**4)
plt.title("y = 2*x^4-x^3-x^2")
def f(x):
    return 2*x**4-x**3-x**2
def h(x):
    return 8*x**3 - 3*x**2 - 2*x
η = 0.05
α = 0.9
v = 0
x = -0.8
iters = 0
X = []
Y = []
while iters<12:
    iters+=1
    X.append(x)
    Y.append(f(x))
    v = α*v - η*h(x)
    x = x + v
    print(iters,x)
plt.plot(X,Y)
plt.show()

import tensorflow as tf
tf.compat.v1.train.MomentumOptimizer(self, learning_rate, momentum,use_locking=False, name="Momentum", use_nesterov=False)

7.Nesterov accelerated gradient (NAG)

由上面Momentum公式，可以畫出Momentum第一和第二次參數更新的方向：

Momentum的想法很簡單，就是多更新一部分上一次迭代的更新量，來平滑這一次迭代的梯度。從物理的角度上解釋，就像是一個小球滾落的時候會受到自身歷史動量的影響，所以才叫動量（Momentum）算法。這樣做直接的效果就是使得梯度下降的的時候轉彎掉頭的幅度不那么大了，於是就能夠更加平穩、快速地沖向局部最小點。

然后NAG就對Momentum說：“既然我都知道我這一次一定會走αv_t-1的量，那么我何必還用現在這個位置的梯度呢？我直接先走到αv_t-1之后的地方，然后再根據那里的梯度再前進一下，豈不美哉？”所以就有了下面的公式：

跟上面Momentum公式的唯一區別在於，梯度不是根據當前參數位置W_t，而是根據先走了本來計划要走的一步后，達到的參數位置W_t - αv_t-1計算出來的。

 

動量解決SGD的兩個問題：（1）SGD引入的噪聲（2）Hessian矩陣病態（SGD收斂過程的梯度相比正常來回擺動幅度較大，NAG是Momentum變種，Nesterov accelerated gradient（NAG，涅斯捷羅夫梯度加速）不僅增加了動量項，並且在計算參數的梯度時，在損失函數中減去了動量項，即計算∇θJ(θ−γνt−1)，這種方式預估了下一次參數所在的位置。 NAG的計算在模型參數施加當前速度之后，可以理解為在Momentum 中引入了一個校正因子。在Momentum中，小球會盲目的跟從下坡的梯度，易發生錯誤。因此，需要提前知道下降的方向，同時，在快到目標點時速度會有所下降，以不至於超出。    可以表示小球下一個大概的位置，從而知道下一個位置的梯度，然后使用當前位置來更新參數。

優點： 

這種基於預測的更新方法，使我們避免過快地前進，並提高了算法地響應能力，NAG 可以使 RNN 在很多任務上有更好的表現。大大改進了 RNN 在一些任務上的表現【為什么對RNN好呢，不懂啊】 

沒有對比就沒有傷害，NAG方法收斂速度明顯加快。波動也小了很多。實際上NAG方法用到了二階信息，所以才會有這么好的結果。先按照原來的梯度走一步的時候已經求了一次梯度，后面再修正的時候又求了一次梯度，所以是二階信息。 

為了從另一個角度更加深入地理解這個算法，我們可以對NAG原來的更新公式進行變換，得到這樣的等效形式:

NAG的原始形式到等效形式的推導。由

其中，和分別是這一次和上一次的更新方向，表示目標函數在處的梯度，超參數是對上一次更新方向的衰減權重，所以一般是0到1之間，是學習率。

可得

記:

上式代入上上式，就得到了NAG等效形式的第二個式子：

對展開可得:

於是我們可以寫出的形式，然后用減去消去后面的無窮多項，就得到了NAG等效形式的第一個式子：

最終我們就得到了NAG的等效形式：

這個NAG的等效形式與Momentum的區別在於，本次更新方向多加了一個，它的直觀含義就很明顯了：如果這次的梯度比上次的梯度變大了，那么有理由相信它會繼續變大下去，那我就把預計要增大的部分提前加進來；如果相比上次變小了，也是類似的情況。這樣的解釋聽起來好像和原本的解釋一樣玄，但是讀者可能已經發現了，這個多加上去的項不就是在近似目標函數的二階導嘛！所以NAG本質上是多考慮了目標函數的二階導信息，怪不得可以加速收斂了！其實所謂“往前看”的說法，在牛頓法這樣的二階方法中也是經常提到的，比喻起來是說“往前看”，數學本質上則是利用了目標函數的二階導信息。

結論：

在原始形式中，Nesterov Accelerated Gradient（NAG）算法相對於Momentum的改進在於，以“向前看”看到的梯度而不是當前位置梯度去更新。經過變換之后的等效形式中，NAG算法相對於Momentum多了一個本次梯度相對上次梯度的變化量，這個變化量本質上是對目標函數二階導的近似。由於利用了二階導的信息，NAG算法才會比Momentum具有更快的收斂速度。

參數設置： 
同Momentum，並把use_nesterov設為True

import tensorflow as tf
tf.compat.v1.train.MomentumOptimizer(learning_rate,momentum,use_nesterov=True)

momentum項和nesterov項都是為了使梯度更新更加靈活，對不同情況有針對性。但是，人工設置一些學習率總還是有些生硬，接下來介紹幾種自適應學習率的方法（AdaGrad，RMSProp，Adam，AdaDelta等。優點：無需人為調節學習率，可以自動調節。缺點；隨着迭代次數增多，學習率越來越小，最終趨於0。）

8.Adagrad

前面的一系列優化算法有一個共同的特點，就是對於每一個參數都用相同的學習率進行更新。但是在實際應用中各個參數的重要性肯定是不一樣的，所以我們對於不同的參數要動態的采取不同的學習率，讓目標函數更快的收斂。 

adagrad方法是將每一個參數的每一次迭代的梯度取平方累加再開方，用基礎學習率除以這個數，來做學習率的動態更新。然后給不同的參數不同的學習率【這樣每一個參數的學習率就與他們的梯度有關系了，那么每一個參數的學習率就不一樣了！也就是所謂的自適應學習率】 

假設N元函數f(x)，針對一個自變量研究Adagrad梯度下降的迭代過程，

 

可以看出，Adagrad算法中有自適應調整梯度的意味（adaptive gradient），學習率需要除以一個東西，這個東西就是前n次迭代過程中偏導數的平方和再加一個常量最后開根號。

舉例：

使用Adagrad算法求y = x²的最小值點

導函數為g(x) = 2x

初始化x⁽⁰⁾ = 4，學習率η=0.25，ε=0.1

from matplotlib import pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
x = np.arange(-4, 4, 0.025)
plt.plot(x,x**2)
plt.title("y = x^2")
def f(x):
    return x**2
def h(x):
    return 2*x
η = 0.25
ε = 0.1
x = 4
iters = 0
sum_square_grad = 0
X = []
Y = []
while iters<12:
    iters+=1
    X.append(x)
    Y.append(f(x))
    sum_square_grad += h(x)**2
    x = x - η/np.sqrt(sum_square_grad+ε)*h(x)
    print(iters,x)
plt.plot(X,Y,"ro")
ax = plt.subplot()
for i in range(len(X)):
    ax.text(X[i], (X[i])**2, "({:.3f},{:.3f})".format(X[i], (X[i])**2), color='red')
plt.show()

公式的另一種寫法：

 其中，⊙是按元素相乘，r為梯度累積變量，r的初始值為0。ε為全局學習率，需要自己設置。δ為小常數，為了防止分母為零。

參數設置： 

• 只需要設置初始學習率，后面學習率會自我調整，越來越小

優點：

• 前期Gt較小的時候， regularizer較大，能夠放大梯度
• 后期Gt較大的時候，regularizer較小，能夠約束梯度
• 學習率將隨着梯度的倒數增長，也就是說較大梯度具有較小的學習率，而較小的梯度具有較大的學習率，可以解決普通的sgd方法中學習率一直不變的問題
• 適合處理稀疏梯度:相當於為每一維參數設定了不同的學習率：壓制常常變化的參數，突出稀缺的更新。能夠更有效地利用少量有意義樣本。在數據分布稀疏的場景上，能更好利用稀疏梯度的信息，比標准的SGD算法更有效地收斂。

缺點： 

• 主要缺陷來自分母項的對梯度平方不斷累積，隨之時間步地增加，分母項越來越大，最終導致學習率收縮到太小無法進行有效更新。如果初始梯度很大的話，會導致整個訓練過程的學習率一直很小，從而導致學習時間變長。

import tensorflow as tf
tf.compat.v1.train.AdagradOptimizer()
tf.compat.v1.train.ProximalAdagradOptimizer() #用來優化目標函數中存在有不可部分的或是加了L1或L2正則的。
tf.compat.v1.train.AdagradDAOptimizer() 

注意：AdagradDA全稱是Adagrad雙重平均算法（Adagrad Dual Averaging algorithm），這個算法我實在是看不懂了，光是數學原理就要好多頁，這里貼出論文的鏈接http://www.jmlr.org/papers/volume12/duchi11a/duchi11a.pdf

Adagrad的一大優勢時可以避免手動調節學習率，比如設置初始的缺省學習率為0.01，然后就不管它，另其在學習的過程中自己變化。當然它也有缺點，就是它計算時要在分母上計算梯度平方的和，由於所有的梯度的平方必為正數，這樣就造成在訓練的過程中，分母累積的和會越來越大。這樣學習到后來的階段，網絡的更新能力會越來越弱，能學到的更多知識的能力也越來越弱，因為學習率會變得極其小【就會提前停止學習】，為了解決這樣的問題又提出了Adadelta算法。

9.RMSprop

由於AdaGrad算法的機制，導致每個元素的學習率在迭代過程中只能降低或者不變，因此很可能出現早期迭代到不好的極值點之后，由於學習率太小而無法沖出這個極值點導致最后收斂到的解不優，為了解決這一問題，RMSProp是基於AdaGrad算法做了一點小修改，其更新公式為：

 

其中，ϵ是防止分母爆0的常數。

這樣，就有了一個改進版的AdaGrad。

該方法即Tieleman&Hinton的RMSProp，由於RMSProp和AdaDelta是同年出現的，

Matthew D. Zeiler並不知道這種改進的AdaGrad被祖師爺命名了。

RMSProp利用了二階信息做了Gradient優化，在BatchNorm之后，對其需求不是很大。

但是沒有根本實現自適應的學習率，依然需要線性搜索初始學習率，然后對其逐數量級下降。

另外，RMSProp的學習率數值與MomentumSGD差別甚大，需要重新線性搜索初始值。

注意：建議取值為1。
特點：

• 其實RMSprop依然依賴於全局學習率
• RMSprop算是Adagrad的一種發展，和Adadelta的變體，效果趨於二者之間
• 適合處理非平穩目標（也就是與時間有關的）
• 對於RNN效果很好，因為RMSprop的更新只依賴於上一時刻的更新，所以適合。

from matplotlib import pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
x = np.arange(-4, 4, 0.025)
plt.plot(x,x**2)
plt.title("y = x^2")
def f(x):
    return x**2
def h(x):
    return 2*x
g = 1
x = 4
ρ = 0.9
η = 0.01
ε = 10e-10
iters = 0
X = []
Y = []
while iters<12:
    iters+=1
    X.append(x)
    Y.append(f(x))
    g = ρ*g+(1-ρ)*h(x)**2
    x = x - η/np.sqrt(g+ε)*h(x)
    print(iters,x)
plt.plot(X,Y,"ro")
ax = plt.subplot()
for i in range(len(X)):
    ax.text(X[i], (X[i])**2, "({:.3f},{:.3f})".format(X[i], (X[i])**2), color='red')
plt.show()

import tensorflow as tf
tf.compat.v1.train.RMSPropOptimizer(learning_rate)

10.Adadelta

Adagrad會累加之前所有的梯度平方，而Adadelta只累加固定大小的項【其實就是相當於指數滑動平均，只用了前多少步的梯度平方平均值】，並且也不直接存儲這些項，僅僅是近似計算對應的平均值【這也就是指數滑動平均的優點】 

整理得更新公式為：

優點： 

• 不用依賴於全局學習率了 
• 訓練初中期，加速效果不錯，很快 
• 避免參數更新時兩邊單位不統一的問題 

缺點： 

• 訓練后期，反復在局部最小值附近抖動

import tensorflow as tf
tf.compat.v1.train.AdadeltaOptimizer()

11.Adam

Adam = Adaptive + Momentum，顧名思義Adam集成了SGD的一階動量和RMSProp的二階動量。 

 

 就像Adadelta和RMSprop一樣Adam會存儲之前衰減的平方梯度，同時它也會保存之前衰減的梯度。經過一些處理之后再使用類似Adadelta和RMSprop的方式更新參數。

特點：

• 結合了Adagrad善於處理稀疏梯度和RMSprop善於處理非平穩目標的優點
• 對內存需求較小
• 為不同的參數計算不同的自適應學習率
• 也適用於大多非凸優化
• 適用於大數據集和高維空間

from matplotlib import pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
x = np.arange(-4, 4, 0.025)
plt.plot(x,x**2)
plt.title("y = x^2")
def f(x):
    return x**2
def h(x):
    return 2*x
x = 4
m = 0
v = 0
β1 = 0.9
β2 = 0.999
η = 0.001
ε = 10e-8
iters = 0
X = []
Y = []
while iters<12:
    iters+=1
    X.append(x)
    Y.append(f(x))
    m = β1*m + (1-β1)*h(x)
    v = β2*v + (1-β2)*h(x)**2
    m_het = m/(1-β1**iters)
    v_het = v/(1-β2**iters)
    x = x - η/np.sqrt(v_het+ε)*m_het
    print(iters,x)
plt.plot(X,Y,"ro")
ax = plt.subplot()
for i in range(len(X)):
    ax.text(X[i], (X[i])**2, "({:.3f},{:.3f})".format(X[i], (X[i])**2), color='red')
plt.show()

import tensorflow as tf
tf.compat.v1.train.AdamOptimizer()

12.Adamax

13.Nadam

tensorflow的優化器及Keras中的默認參數

import tensorflow as tf
tf.compat.v1.train.GradientDescentOptimizer()
tf.compat.v1.train.AdadeltaOptimizer()
tf.compat.v1.train.AdagradOptimizer()
tf.compat.v1.train.AdagradDAOptimizer()
tf.compat.v1.train.MomentumOptimizer()
tf.compat.v1.train.AdamOptimizer()
tf.compat.v1.train.FtrlOptimizer()
tf.compat.v1.train.ProximalGradientDescentOptimizer()
tf.compat.v1.train.ProximalAdagradOptimizer()
tf.compat.v1.train.RMSPropOptimizer()

tf.keras.optimizers.SGD(lr=0.001,momentum=0.9)
tf.keras.optimizers.Adagrad(lr=0.01,epsilon=1e-8)
tf.keras.optimizers.Adadelta(lr=0.01,rho=0.95,epsilon=1e-8)
tf.keras.optimizers.RMSprop(lr=0.001,rho=0.9,epsilon=1e-8)
tf.keras.optimizers.Adam(lr=0.001,beta_1=0.9,beta_2=0.999,epsilon=1e-8)

tensorflow中分布式Tensorflow中同步梯度更新（tf.compat.v1.train.SyncReplicasOptimizer）

在tensorflow的訓練中，分布式可以大大的加快模型訓練速度，但是分布式怎么分配和參數設定，都和SyncReplicasOptimizer這個函數有很大關系。

首先這個SyncReplicasOptimizer函數是專門出來分布式深度學習中的同步梯度下降的(要用異步梯度的可以直接無視)

def __init__(self,
opt: Any,
replicas_to_aggregate: Any,
total_num_replicas: Any = None,
variable_averages: Any = None,
variables_to_average: Any = None,
use_locking: bool = False,
name: str = "sync_replicas") -> None
Construct a sync_replicas optimizer.
Params:
opt – The actual optimizer that will be used to compute and apply the gradients. Must be one of the Optimizer classes.
replicas_to_aggregate – number of replicas to aggregate for each variable update.
total_num_replicas – Total number of tasks/workers/replicas, could be different from replicas_to_aggregate. If total_num_replicas > replicas_to_aggregate: it is backup_replicas + replicas_to_aggregate. If total_num_replicas replicas_to_aggregate: Replicas compute multiple batches per update to variables.
variable_averages – Optional `ExponentialMovingAverage` object, used to maintain moving averages for the variables passed in `variables_to_average`.
variables_to_average – a list of variables that need to be averaged. Only needed if variable_averages is passed in.
use_locking – If True use locks for update operation.
name – string. Optional name of the returned operation.

 opt：將用於計算和應用梯度的實際優化器。必須是Optimizer類之一。
 variable_averages：可選的`ExponentialMovingAverage`對象，用於保持傳入的變量的移動平均值
variables_to_average：需要平均的變量列表。只要如果傳入variable_averages則需要。
use_locking：如果True使用鎖定進行更新操作。
name：string。返回操作的可選名稱。

除開opt之外，其他可以先不指定，我們這里重點看看replicas_to_aggregate和total_num_replicas這兩個參數，這是非常重要的。

replicas_to_aggregate：是指我們所需要通過分布式的worker計算完成之后所需要的副本數，這里可以認為就是梯度的個數

total_num_replicas：是指我們指定出入計算樣本的數量（幾個機器就分配幾個batch_size） 

#同步模式計算更新梯度
rep_op = tf.train.SyncReplicasOptimizer(optimizer,
                                        replicas_to_aggregate=len(worker_hosts),
                                        total_num_replicas=len(worker_hosts),
                                        use_locking=True)

我們可以看出，replicas_to_aggregate和 total_num_replicas都是等於集群中worker的數量的，但是源碼中對於這兩個參數又添加了額外的解釋，(replicas_to_aggregate以下簡稱副本數，total_num_replicas以下簡稱計算數)

replicas_to_aggregate: number of replicas to aggregate for each variable
update.
total_num_replicas: Total number of tasks/workers/replicas, could be
different from replicas_to_aggregate.
If total_num_replicas > replicas_to_aggregate: it is backup_replicas +
replicas_to_aggregate.
If total_num_replicas < replicas_to_aggregate: Replicas compute
multiple batches per update to variables.

 什么意思呢？意思就是說，副本數可以和計算數可以是不相等的。就好比有三個人交了三張試卷一樣，三個人就是計算數，三張試卷就是副本數。

 源碼中說，If total_num_replicas > replicas_to_aggregate: it is backup_replicas + replicas_to_aggregate.

這是指如果計算數大於副本數，那么多的那個就是backup worker，即備份工作節點。

比如：四個人交了三張試卷，沒交試卷的我們就認為是備份工作節點   

那這又是什么玩意兒呢，圖啥呢？

在同步梯度下降中，所有worker都要等大家一起計算完成之后才能進行計算，因此會出現等待機制。就像四個人做試卷，但是每個人的速度不一樣，但是要同時交卷，所以我們要等待交卷最慢的那個人算完之后才能算完，這樣很浪費算的很快的那個人的能力，所以就有了backup的提出，4個人交卷我們只要三個人的試卷，最慢的那個我們不要了，直接清除掉，開始下一輪的考試~

因此，加入backup worker的方法可以加速分布式同步梯度下法的訓練模型，那么backup取幾個呢？Revisiting Distributed Synchronous SGD這篇論文中認為是20：1的比例。

還有一種情況，那就是計算數大於副本數。當計算數小於副本數，那么他會多計算其他worker的樣本，就好比三個人(給三人一人發一個試卷，分配的是沒多余的)，但是要交四份試卷，那就先各算各的，先計算完成的再從已知的三份試卷中拿一個計算，然后再等着這三個人一共上交了四份試卷之后，這樣就ok了。有什么用呢？我也不清楚，我們看下一步

當我們說這計算數和副本數的時候總是和機器的物理數量有關系（gpu數量），那么其實這兩個都可以不等於機器的物理數量。例如：

gpu數量：3           那么len（worker）就是3

計算數：10           那么這10就指的是一次給分10個batch_size，然后讓3個worker去算

副本數：8             那么這是個batch_size算完之后我只要8個梯度，然后平均，參數服務器的什么后續操作

import numpy as np
import tensorflow as tf

F = tf.flags.FLAGS
# Define parameters
FLAGS = tf.flags.FLAGS
tf.flags.DEFINE_float('learning_rate', 0.00003, 'Initial learning rate.')
tf.flags.DEFINE_integer('steps_to_validate', 1000,
                        'Steps to validate and print loss')

# For distributed
tf.flags.DEFINE_string("ps_hosts", "",
                       "Comma-separated list of hostname:port pairs")
tf.flags.DEFINE_string("worker_hosts", "",
                       "Comma-separated list of hostname:port pairs")
tf.flags.DEFINE_string("job_name", "", "One of 'ps', 'worker'")
tf.flags.DEFINE_integer("task_index", 0, "Index of task within the job")
tf.flags.DEFINE_integer("issync", 0, "是否采用分布式的同步模式，1表示同步模式，0表示異步模式")

# Hyperparameters
learning_rate = FLAGS.learning_rate
steps_to_validate = FLAGS.steps_to_validate

def main(_):
    ps_hosts = FLAGS.ps_hosts.split(",")
    worker_hosts = FLAGS.worker_hosts.split(",")
    cluster = tf.train.ClusterSpec({"ps": ps_hosts, "worker": worker_hosts})
    server = tf.train.Server(cluster, job_name=FLAGS.job_name, task_index=FLAGS.task_index)

    issync = FLAGS.issync
    if FLAGS.job_name == "ps":
        server.join()
    elif FLAGS.job_name == "worker":
        with tf.device(tf.train.replica_device_setter(
                worker_device="/job:worker/task:%d" % FLAGS.task_index,
                cluster=cluster)):
            global_step = tf.Variable(0, name='global_step', trainable=False)

            input = tf.placeholder("float")
            label = tf.placeholder("float")

            weight = tf.get_variable("weight", [1], tf.float32, initializer=tf.random_normal_initializer())
            biase = tf.get_variable("biase", [1], tf.float32, initializer=tf.random_normal_initializer())
            pred = tf.matmul(input, weight) + biase

            loss_value = loss(label, pred)
            optimizer = tf.train.GradientDescentOptimizer(learning_rate)

            grads_and_vars = optimizer.compute_gradients(loss_value)
            if issync == 1:
                # 同步模式計算更新梯度
                rep_op = tf.train.SyncReplicasOptimizer(optimizer,
                                                        replicas_to_aggregate=len(
                                                            worker_hosts),
                                                        replica_id=FLAGS.task_index,
                                                        total_num_replicas=len(
                                                            worker_hosts),
                                                        use_locking=True)
                train_op = rep_op.apply_gradients(grads_and_vars,
                                                  global_step=global_step)
                init_token_op = rep_op.get_init_tokens_op()
                chief_queue_runner = rep_op.get_chief_queue_runner()
            else:
                # 異步模式計算更新梯度
                train_op = optimizer.apply_gradients(grads_and_vars,
                                                     global_step=global_step)

            init_op = tf.initialize_all_variables()

            saver = tf.compat.v1.train.Saver()
            tf.summary.scalar('cost', loss_value)
            summary_op = tf.summary.merge_all()

        sv = tf.train.Supervisor(is_chief=(FLAGS.task_index == 0),
                                 logdir="./checkpoint/",
                                 init_op=init_op,
                                 summary_op=None,
                                 saver=saver,
                                 global_step=global_step,
                                 save_model_secs=60)

        with sv.prepare_or_wait_for_session(server.target) as sess:
            # 如果是同步模式
            if FLAGS.task_index == 0 and issync == 1:
                sv.start_queue_runners(sess, [chief_queue_runner])
                sess.run(init_token_op)
            step = 0
            while step < 1000000:
                train_x = np.random.randn(1)
                train_y = 2 * train_x + np.random.randn(1) * 0.33 + 10
                _, loss_v, step = sess.run([train_op, loss_value, global_step],
                                           feed_dict={input: train_x, label: train_y})
                if step % steps_to_validate == 0:
                    w, b = sess.run([weight, biase])
                    print("step: %d, weight: %f, biase: %f, loss: %f" % (step, w, b, loss_v))

        sv.stop()


def loss(label, pred):
    return tf.square(label - pred)


if __name__ == "__main__":
    tf.app.run()

14.牛頓法——二階優化方法

1）牛頓法（Newton's method）

牛頓法是一種在實數域和復數域上近似求解方程的方法。方法使用函數f (x)的泰勒級數的前面幾項來尋找方程f (x) = 0的根。牛頓法最大的特點就在於它的收斂速度很快。

具體步驟：

首先，選擇一個接近函數 f (x)零點的 x₀，計算相應的 f (x₀) 和切線斜率f  ' (x₀)（這里f ' 表示函數 f  的導數）。然后我們計算穿過點(x_0,  f  (x₀)) 並且斜率為f '(x₀)的直線和 x 軸的交點的x坐標，也就是求如下方程的解：

我們將新求得的點的 x 坐標命名為x₁，通常x₁會比x₀更接近方程f  (x) = 0的解。因此我們現在可以利用x₁開始下一輪迭代。迭代公式可化簡為如下所示：

已經證明，如果f  ' 是連續的，並且待求的零點x是孤立的，那么在零點x周圍存在一個區域，只要初始值x₀位於這個鄰近區域內，那么牛頓法必定收斂。 並且，如果f  ' (x)不為0, 那么牛頓法將具有平方收斂的性能. 粗略的說，這意味着每迭代一次，牛頓法結果的有效數字將增加一倍。下圖為一個牛頓法執行過程的例子。

　　由於牛頓法是基於當前位置的切線來確定下一次的位置，所以牛頓法又被很形象地稱為是"切線法"。牛頓法的搜索路徑（二維情況）如下圖所示：

　　牛頓法搜索動態示例圖：

 

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

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

 

第二步，把 x 看做自變量， 所有帶有 x^k 的項看做常量，令一階導數為 0 ，即可求近似函數的最小值：

 

即：

 

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

牛頓法是二階算法，因為該算法使用了海森矩陣（Hessian matrix）求權重的二階偏導數。牛頓法的目標就是采用損失函數的二階偏導數尋找更好的訓練方向。

因此牛頓下降法是用二次曲面去擬合當前的局部曲面,現在采用如下表示：f(wi) = fi、ᐁf(wi) = gi 和 Hf(wi) = Hi。在w0點使用泰勒級數展開式二次逼近函數f.

H0為函數f在點w0的海森矩陣，對w求導，得到g,通過將g設定為0，我們就可以找到f(w)的最小值，也就得到了以下方程式。

因此，從參數向量w0開始，牛頓法服從以下方式進行迭代：

向量 Hi-1·gi也就是所說的牛頓下降步（Newton's step）。注意，參數的這些變化將朝着極大值而不是極小值逼近，出現這樣的情況是因為海森矩陣非正定。

因此在不能保證矩陣正定的情況下，損失函數並不能保證在每一次迭代中都是減少的。為了防止上述問題，牛頓法的方程式通常可以修改為：

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

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

import numpy as np
def Jacobian(x):
    return np.array([x[0], 0.4 * x[1], 1.2 * x[2]])

def Hessian(x):
    return np.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 = np.dot(np.linalg.inv(Hessian(x)), Jacobian(x))
        xOld = x
        x = x + alpha * p
        Delta = sum((x - xOld) ** 2)
        i += 1
    print(x)
x0 = np.array([-2, 2, -2])
Newton(x0)

關於牛頓法和梯度下降法的效率對比：

從本質上去看，牛頓法是二階收斂，梯度下降是一階收斂，所以牛頓法就更快。如果更通俗地說的話，比如你想找一條最短的路徑走到一個盆地的最底部，梯度下降法每次只從你當前所處位置選一個坡度最大的方向走一步，牛頓法在選擇方向時，不僅會考慮坡度是否夠大，還會考慮你走了一步之后，坡度是否會變得更大。所以，可以說牛頓法比梯度下降法看得更遠一點，能更快地走到最底部。（牛頓法目光更加長遠，所以少走彎路；相對而言，梯度下降法只考慮了局部的最優，沒有全局思想。）

根據wiki上的解釋，從幾何上說，牛頓法就是用一個二次曲面去擬合你當前所處位置的局部曲面，而梯度下降法是用一個平面去擬合當前的局部曲面，通常情況下，二次曲面的擬合會比平面更好，所以牛頓法選擇的下降路徑會更符合真實的最優下降路徑。

 

注：紅色的牛頓法的迭代路徑，綠色的是梯度下降法的迭代路徑。

牛頓法的優缺點總結：

• 優點：二階收斂，收斂速度快（迭代次數少）；
• 缺點：牛頓法是一種迭代算法，每一步都需要求解目標函數的Hessian矩陣的逆矩陣，計算比較復雜。

15.擬牛頓法（Quasi-Newton method）

擬牛頓法是求解非線性優化問題最有效的方法之一，於20世紀50年代由美國Argonne國家實驗室的物理學家W.C.Davidon所提出來。Davidon設計的這種算法在當時看來是非線性優化領域最具創造性的發明之一。不久R. Fletcher和M. J. D. Powell證實了這種新的算法遠比其他方法快速和可靠，使得非線性優化這門學科在一夜之間突飛猛進。

擬牛頓法的本質思想是改善牛頓法每次需要求解復雜的Hessian矩陣的逆矩陣的缺陷，它使用正定矩陣來近似Hessian矩陣的逆，從而簡化了運算的復雜度。擬牛頓法和最速下降法一樣只要求每一步迭代時知道目標函數的梯度。通過測量梯度的變化，構造一個目標函數的模型使之足以產生超線性收斂性。這類方法大大優於最速下降法，尤其對於困難的問題。另外，因為擬牛頓法不需要二階導數的信息，所以有時比牛頓法更為有效。如今，優化軟件中包含了大量的擬牛頓算法用來解決無約束，約束，和大規模的優化問題。

具體步驟：

擬牛頓法的基本思想如下。首先構造目標函數在當前迭代x_k的二次模型：

這里B _k是一個對稱正定矩陣，於是我們取這個二次模型的最優解作為搜索方向，並且得到新的迭代點：

其中我們要求步長a_k 滿足Wolfe條件。這樣的迭代與牛頓法類似，區別就在於用近似的Hesse矩陣B_k   代替真實的Hesse矩陣。所以擬牛頓法最關鍵的地方就是每一步迭代中矩陣B_k的更新。現在假設得到一個新的迭代x_k+1，並得到一個新的二次模型：

我們盡可能地利用上一步的信息來選取B _k。具體地，我們要求  

從而得到

這個公式被稱為割線方程。 常用的擬牛頓法有DFP算法和BFGS算法。

因為需要很多的操作求解海森矩陣的值還有計算矩陣的逆，應用牛頓法的計算量是十分巨大的。因此有一種稱之為擬牛頓法（quasi-Newton）或

變量矩陣法來解決這樣的缺點。這些方法並不是直接計算海森矩陣然后求其矩陣的逆，擬牛頓法是在每次迭代的時候計算一個矩陣，其逼近海森矩陣的逆。

最重要的是，該逼近值只是使用損失函數的一階偏導來計算。

海森矩陣由損失函數的二階偏導組成，擬牛頓法背后的思想主要是僅使用損失函數的一階偏導數，通過另一矩陣G逼近海森矩陣的逆。擬牛頓法的公式可以表示為：

學習速率η可以設定為固定常數，也可以通過單變量函數優化得到。其中矩陣G逼近海森矩陣的逆，且有不同的方式進行逼近。通常，最常用的兩種方式是Davidon–Fletcher–Powell formula (DFP) 和 the Broyden–Fletcher–Goldfarb–Shanno formula (BFGS)。

擬牛頓法的訓練過程流程圖如下所示，從圖中我們可以看出來模型是通過第一次計算擬牛頓訓練方向而優化參數的，然后再尋找適當的學習速率。
擬牛頓法適用於絕大多數案例中：它比梯度下降和共軛梯度法收斂更快，並且也不需要確切地計算海森矩陣及其逆矩陣。

16.Levenberg–Marquardt Algorithm

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

Levenberg-Marquardt 算法，也稱為衰減最小二乘法（damped least-squares method），該算法的損失函數采用平方誤差和的形式。該算法的執行也不需要

計算具體的海森矩陣，它僅僅只是使用梯度向量和雅克比矩陣（Jacobian matrix）。

這里解釋一下雅克比矩陣（Jacobian matrix）：Jacobi 矩陣實際上是向量值函數的梯度矩陣，假設F:Rn→Rm 是一個從n維歐氏空間轉換到m維歐氏空間的函數。

這個函數由m個實函數組成:

這些函數的偏導數(如果存在)可以組成一個m行n列的矩陣(m by n)，這就是所謂的雅可比矩陣：

該算法的損失函數如下方程式所示為平方誤差和的形式：

在上式中，m是數據集樣本的數量。
我們可以定義損失函數的雅可比矩陣以誤差對參數的偏導數為元素，公式如下：

其中m是數據集樣本的數量，n是神經網絡的參數向量。那么雅可比矩陣是m*n階矩陣，損失函數的梯度向量就可以按如下計算出來：

E在這里是所有誤差項的向量。
最終，我們可以利用一下表達式逼近海森矩陣：

其中λ為衰減因子，它確保了海森矩陣的證定性（positiveness），I是單位矩陣。下面的表達式定義了Levenberg-Marquardt 算法中參數的更新和優化過程：

當衰減參數λ為0時，Levenberg-Marquardt 算法就是使用海塞矩陣逼近值的牛頓法。而當λ很大時，該算法就近似於采用很小學習速率的梯度下降法。如果進行迭代導致了

損失函數上升，衰減因子λ就會下降，從而Levenberg-Marquardt 算法更接近於牛頓法。該過程經常用於加速收斂到極小值點。

第一步就是計算損失函數、梯度和海森矩陣逼近值，隨后再在每次迭代降低損失中調整衰減參數。

Levenberg-Marquardt 算法是為平方誤差和函數所定制的。這就讓使用這種誤差度量的神經網絡訓練地十分迅速。但是其缺點是：

• 1. 它不能用於平方根誤差或交叉熵誤差（cross entropy error）等函數。
• 2. 該算法與正則項不兼容。
• 3. 對於大型數據集或神經網絡，雅可比矩陣變得十分巨大，因此也需要大量的內存。
• 在大型數據集或神經網絡中並不推薦采用Levenberg-Marquardt 算法。

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


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

        pnew = np.zeros(np.shape(p))
        nits2 = 0
        while (p != pnew).all() and nits2 < maxits:
            nits2 += 1
            dp, resid, rank, s = np.linalg.lstsq(H, grad)
            pnew = p - dp
            fpnew, rnew, gradnew, Jnew = function(pnew)
            enew = np.sum(np.dot(np.transpose(rnew), rnew))
            rho = np.linalg.norm(np.dot(np.transpose(r), r) - np.dot(np.transpose(rnew), rnew))
            rho /= np.linalg.norm(np.dot(np.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, np.linalg.norm(grad), nu)
        
p0 = np.array([-1.92, 2])
lm(p0)

17.共軛梯度法（Conjugate gradient）

共軛方向法是介於最速下降法和Newton法之間的一種方法。克服了最速下降法的鋸齒現象，從而提高了收斂速度；同時，共軛方向法的迭代公式比較簡單，不必求目標函數的Hesse矩陣，比Newton法減少了計算量和存儲量。是一種比較實用而且有效的方法。

共軛梯度法是共軛方向法的一種，意思是搜索方向都互相共軛

共軛向量及其性質

定義1：(共軛方向) Q是n×n 對稱正定矩陣，若 nn 維向量空間中的非零向量 p0，p1，⋅⋅⋅，pm−1p0，p1，···，pm−1 滿足

則稱p0，p1，⋅⋅⋅，pm−1p0，p1，···，pm−1 是Q共軛向量或稱向量p0，p1，⋅⋅⋅，pm−1p0，p1，···，pm−1是Q共軛的(簡稱共軛)，稱p0，p1，⋅⋅⋅，pm−1p0，p1，···，pm−1的方向是Q共軛方向。
當Q=I(單位矩陣)Q=I(單位矩陣)時，公式(1)變為

即向量p0，p1，⋅⋅⋅，pm−1互相正交。由此可見，正交是共軛的一種特殊情況，共軛是正交的推廣。
定理1： 若非零向量p0，p1，⋅⋅⋅，pm−1是Q共軛的，則他們是線性無關的。
推論1： 在n維的向量空間中，非零的共軛向量個數不超過n。
定義2： 設p0，p1，⋅⋅⋅，pm−1是Rn中線性無關向量，x0∈Rⁿ，Pⁱ∈Rⁿ那么由形式為

z的向量構成的集合為由點x0和向量p0，p1，⋅⋅⋅，pm−1所生成的線性流行。記為L[x0;p0,p1,···,pm−1]。

基本思想

在考慮普通函數之前，我們首先用2元正定二次函數進行講解。首先考慮如下的正定二次函數

要求的目標函數f(x)的最優值，根據最速下降法的思想，我們首先選定一個初始點x0x0，然后沿着該點的最速下降方向p0=−∇f(x0)做直線搜索，得到點x1，由最速下降法的性質(可以想一下最速下降法中的鋸齒現象)可知

在第二次迭代過程中我們不適用−∇f(x1)作為這次迭代的搜索方向，我們想直在第二次迭代之后能直接到達最優點x^*（x^*表示是當前迭代找到的最優點），那么這次的迭代方向p1應該滿足什么條件呢？
首先根據最速下降法的迭代公式我們有

其中t₁是最優步長因子,顯然在未到達最優點x^*之前，t₁是不等於0的。對目標函數求梯度，有梯度函數

對於最優點x^*，我們有

即

在上式兩邊同乘以P₀^T，由於公式(4)，並且t₁≠0我們可以得到

 

由公式(7)我們知道，p0p0和p1p1為QQ的共軛向量。 
現在我們假設:

 

 

 

綜上所述，對於n元正定二次函數，我們可以從任意點出發，然后沿着這n個共軛方向最多做n次直線搜索(原因將在下一小節進行講解)，就可以求的目標函數的最優點。

優化模型建立

雖然梯度下降法的每一步都是朝着局部最優的方向前進的，但是它在不同的迭代輪數中會選擇非常近似的方向，說明這個方向的誤差並沒通過一次更新方向和步長更新完，在這個方向上還存在誤差，因此參數更新的軌跡是鋸齒狀。共軛梯度法的思想是，選擇一個優化方向后，本次選擇的步長能夠將這個方向的誤差更新完，在以后的優化更新過程中不再需要朝這個方向更新了。由於每次將一個方向優化到了極小，后面的優化過程將不再影響之前優化方向上的極小值，所以理論上對N維問題求極小只用對N個方向都求出極小就行了。為了不影響之前優化方向上的更新量，需要每次優化方向共軛正交。假定每一步的優化方向用p_k表示，可得共軛正交

由此可得，每一步優化后，當前的誤差和剛才的優化方向共軛正交。

若為N維空間優化問題，則每次優化方向可以組成這個空間中的一組基底。P={p₁,p2,…,p_n} 。

算法只需要解決兩個問題：

1. 優化方向
2. 優化步長

優化方向確定

優化步長的選取

三個推論

最終簡化

算法在三中基本推導完畢，但是在工程應用中如果每次進行p_k的正交化需要對之前所有的優化向量求解β，現簡化如下：

最終的推導結束。整理為如下的偽代碼（用於正定二次函數的共軛梯度法）

對於正定二次函數(3)，從任意點出發，順次沿nn個QQ共軛方向做直線搜索，最多經過nn次就可以到達極小點。 
下面是對於正定二次函數的共軛方向法的算法描述。

import numpy as np
np.random.seed(12)
A = np.random.randint(8,18,size=(2,2))
A = A+A.T #把A變成對稱矩陣
B = np.random.randint(-5,8,size=(2,1))

def func(X):
    return 0.5*np.dot(X.T,np.dot(A,X)) - np.dot(B.T,X)

def d_fun(X):
    return np.dot(A,X) - B

#最速下降法
def gd(x,func,d_f,lr=0.005,e=0.0001):
    init_x = x
    init_y = func(init_x)
    ls = 1
    x_point = init_x
    y_point = init_y

    while ls > e:
        init_x -= lr*d_f(init_x)
        new_y = func(init_x)
        ls = np.abs(new_y-init_y)[0,0]
        init_y = new_y
        x_point = np.hstack((x_point,init_x))
        y_point = np.hstack((y_point,init_y))
    return x_point.T,y_point.T

import copy
#共軛下降法
def cad(init_x,d_f):
    print(init_x)
    x_0 = copy.deepcopy(init_x)
    r_0 = -d_f(x_0)
    p_0 = r_0
    r_k = r_0
    p_k = p_0
    x_k = x_0
    k = 1
    n = len(x_0)
    print(n)
    x_point = x_k
    while True:
        alpha_k = np.dot(r_k.T,r_k)[0,0]/np.dot(p_k.T,np.dot(A,p_k))[0,0]
        x_k += alpha_k*p_k
        r_k2 = r_k - alpha_k*np.dot(A,p_k)
        x_point = np.hstack((x_point,x_k))
        if k >= n:
            break
        beta_k = np.dot(r_k2.T,r_k2)[0,0]/np.dot(r_k.T,r_k)[0,0]
        p_k = r_k2 + beta_k*p_k
        r_k = r_k2
        k += 1
    x_point[:,0] = init_x.reshape(-1)
    return x_point.T

x_p = cad(np.array([[0],[0]],dtype=np.float),d_fun)
x_list,y_list = gd(np.array([[0],[0]],dtype=np.float),func,d_fun)
x1_max, x2_max = np.max(x_list,axis=0)
x1_min, x2_min = np.min(x_list,axis=0)
x1 = np.linspace(x1_min-1,x1_max+1,100)
x2 = np.linspace(x2_min-1, x2_max+1,100)
X1,X2 = np.meshgrid(x1,x2)
X = np.array(list(zip(X1.flatten(),X2.flatten())))
Y = np.array(list(map(lambda x:func(x.reshape(2,1)),X))).reshape(*X1.shape)
contour = plt.contour(X1,X2,Y,levels=list(np.linspace(np.min(y_list)+0.2,np.max(y_list),5)),colors='k')
plt.clabel(contour)
plt.scatter(x_p[:,0],x_p[:,1],c='g')
plt.plot(x_p[:,0],x_p[:,1],c='g',label='cad')
plt.plot(x_list[:,0],x_list[:,1],'r-',label='gd')
plt.legend()
plt.show()

在算法描述的第4步中，只是說要求得共軛方向，並沒有說明如何求共軛方向，不同求共軛方向的方法對應了不同的共軛方向法。所以可以按照上述方法，依次構造出共軛向量，最多經過n次迭代就能找到最優點x^*。 
共軛梯度法的算法描述如下。

用於非二次函數的共軛梯度法

用於一般函數的Fletcher-Reeves共軛梯度法描述如下。

總結

共軛梯度法是一種典型的共軛方向法，它的搜索方向是負梯度方向和上一次搜索方向的一個組合

關於β_k-1，有兩種常用的計算公式

（PRP公式）

（PRP公式）

預處理共軛法改善了G的條件數，使算法的收斂速度加快

預處理的方法是尋找一個非奇異矩陣C，使得C^-TGC^-1的條件數小於G的條件數

再開始共軛梯度法是滿足某一條件后重新使用最速下降方向作為搜索方向，這個條件包括迭代n步，或共軛梯度方向是上升方向

FR共軛梯度法：

# coding=utf-8
# FR-CG共軛梯度法
from linear_search.Function import *
from linear_search.wolfe import *
import numpy as np


def conjugate(_f, _x):
    fun = Function(_f)
    x = array(_x)
    d = -fun.grad(x)
    while fun.norm(x) > 0.01:
        alpha = wolfe(_f, x, d)
        g = np.mat(fun.grad(x))
        beta = 1 / np.dot(g, g.T)
        x = x + alpha * d
        g = np.mat(fun.grad(x))
        beta = beta * np.dot(g, g.T)
        d = array(-g + beta * d)[0]
    return x

 再開始共軛梯度法：

# coding=utf-8
# 再開始共軛梯度法
from linear_search.Function import *
from linear_search.wolfe import *
import numpy as np


def restart_conjugate(_f, _x, n):
    fun = Function(_f)
    x = array(_x)
    while True:
        d = -fun.grad(x)
        k=0
        if(np.linalg.norm(d)<0.01):
            break
        while fun.norm(x) > 0.01:
            g = np.mat(fun.grad(x))

            alpha = wolfe(_f, x, d)
            x = x + alpha * d
            k=k+1

            g1= np.mat(fun.grad(x))

            if np.dot(g1,g.T)/np.dot(g1,g1.T)>0.1 or k>=n:
                if np.linalg.norm(g1)<0.01:
                    return x
                break
            else:
                beta = np.dot(g1, g1.T) / np.dot(g, g.T)
                d = array(-g + beta * d)[0]
                if np.dot(mat(d),g1.T)>0:
                    break
    return x

其他的文檔介紹

共軛梯度法是介於最速下降法與牛頓法之間的一個方法，它僅需利用一階導數信息，但克服了最速下降法收斂慢的缺點，又避免了牛頓法需要存儲和計算Hesse矩陣並求逆的缺點，該算法希望能加速梯度下降的收斂速度，同時避免使用海森矩陣進行求值、存儲和求逆獲得必要的優化信息。共軛梯度法不僅是解決大型線性方程組最有用的方法之一，也是解大型非線性最優化最有效的算法之一。 在各種優化算法中，共軛梯度法是非常重要的一種。其優點是所需存儲量小，具有步收斂性，穩定性高，而且不需要任何外來參數。

在共軛梯度訓練算法中，因為是沿着共軛方向（conjugate directions）執行搜索的，所以通常該算法要比沿着梯度下降方向優化收斂得更迅速。共軛梯度法的訓練方向

是海森矩陣共軛的。我們用d表示訓練方向向量，然后從初始參數向量w0和初始訓練方向向量d0=-g0開始，共軛梯度法所構建的訓練方向序列為：

 

上式中，Y稱之為共軛參數，並且有一些方法計算這個參數。兩種常用的方法源自Fletcher、Reeves 和 Polak 、Ribiere。對於所有的共軛梯度算法，訓練方向

周期性地重置為負梯度方向。參數通過下面的表達式得以更新和優化，通常學習速率η可使用單變量函數優化方法求得。

共軛梯度法的訓練過程流程圖如下表示，從圖中我們可以看出來模型是通過第一次計算共軛梯度訓練方向而優化參數的，然后再尋找適當的學習速率。

共軛梯度法已經證實其在神經網絡中要比梯度下降法有效得多。並且由於共軛梯度法並沒有要求使用海森矩陣，所以在大規模神經網絡中其還可以做到很好的性能。

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

import numpy as np

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

def Hessian(x):
    # return array([[.2,0],[0,1]])
    return np.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 = np.dot(r.transpose(), r)
    beta0 = betaTop

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

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

    while i < iMax and betaTop > epsilon ** 2 * beta0:
        j = 0
        dp = np.dot(p.transpose(), p)
        alpha = (epsilon + 1) ** 2
        # Newton-Raphson iteration
        while j < jMax and alpha ** 2 * dp > epsilon ** 2:
            # Line search
            alpha = np.dot(Jacobian(x).transpose(), p) / (np.dot(p.transpose(), np.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 = np.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 np.dot(r.transpose(), p) <= 0:
            p = r
            k = 0
            print("Restarting")
            
        i += 1
    print(x)
    
x0 = np.array([-2, 2, -2])
CG(x0)

總結

提速可以歸納為以下幾個方面： 

• 使用momentum來保持前進方向(velocity)； 
• 為每一維參數設定不同的學習率：在梯度連續性強的方向上加速前進； 
• 用歷史迭代的平均值歸一化學習率：突出稀有的梯度；

綜上分析，可以得出如下幾個結論：

• AdaGrad、RMSProp、AdaDelta和Adam幾個優化算法，目標函數自變量中每個元素都分別擁有自己的學習率；
• AdaGrad目標函數自變量中各個元素的學習率只能保持下降或者不變，因此當學習率在迭代早期降得較快且當前解依然不佳時，由於后期學習率過小，可能較難找到一個有用的解；
• RMSProp和AdaDelta算法都是解決AdaGrad上述缺點的改進版本，本質思想都是利用最近的時間步的小批量隨機梯度平方項的加權平均來降低學習率，從而使得學習率不是單調遞減的（當最近梯度都較小的時候能夠變大）。不同的是，RMSProp算法還是保留了傳統的學習率超參數，可以顯式指定。而AdaDelta算法沒有顯式的學習率超參數，而是通過ΔxΔx做運算來間接代替學習率；
• Adam算法可以看成是RMSProp算法和動量法的結合。

 

如圖，收斂速度最慢的是梯度下降算法，但該算法同時也只要求最少內存。相反，Levenberg-Marquardt 算法可能是收斂速度最快的，但其同時也要求最多的內存。

比較折中的方法是擬牛頓法。
總而言之，如果我們的神經網絡有數萬參數，為了節約內存，我們可以使用梯度下降或共軛梯度法。如果我們需要訓練多個神經網絡，並且每個神經網絡都只有數百參數、數千樣本，那么我們可以考慮 Levenberg-Marquardt 算法。而其余的情況，擬牛頓法都能很好地應對。

二.在線算法（OGD =》 FOBOS =》 RDA =》 FTRL）

在線學習算法相關背景

在線學習算法強調的是訓練的實時性，面向流式數據，每次訓練不使用全量數據，而是以之前訓練好的參數為基礎，每次利用一個樣本更新一次模型，是屬於增量學習的一部分，從而快速更新模型，提高模型的時效性。

在線學習算法與離線學習關注點的不同

統計學習的先驗假設是數據存在自己一定的分布，我們的目的是尋找與實際分布距離最小的策略來泛化未知的結果。數據由真實模型產生，如果能有無限數據、並在包含有真實模型的空間里求解，也許我們能算出真是模 型。但實際上我們只有有限的有噪音的數據，這又限制我們只能使用相對簡單的模型。所以，理想的算法是能夠用不多的數據來得到一個不錯的模型。

離線學習算法強調的是數據中不斷訓練已有的模型，不斷迭代達到真實的數據分布，在訓練過程中，所有數據都是可見的，目標是一定的（就是最終的那個真實分布），其中可以采用不同的更新策略和采樣策略，所以有了批量梯度下降和隨機梯度下降等算法。

但是，在線學習算法的關注點是不同的，在線學習算法的限制條件是只能看到過去的數據和當前的數據，看不到未來的數據，所以我們訓練的策略可能是任意的數據分布，通過不斷的訓練可能顛覆之前的訓練結果，所以在線學習甚至是為了破壞我們之前的策略精心設計的。

在線學習關注點在於，追求對所知道的所有知識所能設計最優的策略，那么同這個最優的策略的差距成為后悔（regret）:后悔沒有從一開始就選擇這個策略，當然，我們希望的是，隨着時間的增加，這個差異會不斷的變小。因為我們不對數據進行任何假設，所以策略是否完美並不是我們關心的（比如回答所有問題），我們追求的是，沒有后悔（no-regret）

batch模式和delta模式

梯度下降可以分成兩種模式，batch模式和delta模式。batch模式的時效性比delta模式要低一些。分析一下batch模式，比如昨天及昨天的數據訓練成了模型M，那么今天的每一條訓練數據在訓練過程中都會更新一次模型M，從而生成今天的模型M1。

batch學習可以認為是離線學習算法，強調的是每次訓練都需要使用全量的樣本，因而可能會面臨數據量過大的問題。一般進行多輪迭代來向最優解靠近。online learning沒有多輪的概念，如果數據量不夠或訓練數據不夠充分，通過copy多份同樣的訓練數據來模擬batch learning的多輪訓練也是有效的方法。

delta模式可以認為是在線學習算法，沒有多輪的概念，如果數據量不夠或訓練數據不夠充分，通過copy多份同樣的訓練數據來模擬batch learning的多輪訓練也是有效的方法。所以，OGD和SGD都屬於在線學習算法，因為每次更新模型只用一個樣本。SGD則每次只針對一個觀測到的樣本進行更新。通常情況下，SGD能夠比GD“更快”地令 逼近最優值。當樣本數 特別大的時候，SGD的優勢更加明顯，並且由於SGD針對觀測到的“一條”樣本更新 ，很適合進行增量計算，實現梯度下降的Online模式（OGD, OnlineGradient Descent）。

在線凸優化

現在做在線學習和CTR常常會用到邏輯回歸（ Logistic Regression），而傳統的批量（batch）算法無法有效地處理超大規模的數據集和在線數據流，google先后三年時間（2010年-2013年）從理論研究到實際工程化實現的FTRL（Follow-the-regularized-Leader）算法，在處理諸如邏輯回歸之類的帶非光滑正則化項（例如1范數，做模型復雜度控制和稀疏化）的凸優化問題上性能非常出色，據聞國內各大互聯網公司都第一時間應用到了實際產品中，我們的系統也使用了該算法。這里對FTRL相關發展背景和工程實現的一些指導點做一些介紹，凸優化的理論細節不做詳細介紹，感興趣可以去查閱相應paper，相關paper列表會在文后附上。機器學習並非本人在校時的專業方向，不過在校期間積累的基礎不算太差，而且很多東西也是相通的，鑽研一下基本意思都還能搞明白。當然，有不准確的地方歡迎大家討論指正。

本文主要會分三個部分介紹，如果對理論產生背景不感興趣的話，可以直接看第3部分的工程實現（這一部分google13年那篇工程化的paper介紹得很詳細）：

1. 相關背景：包括通用性的問題描述、批量算法、傳統在線學習算法等
2. 簡單介紹與FTRL關系比較密切的Truncated Gradient、FOBOS以及RDA（Regularized Dual Averaging）等算法
3. FTRL理論公式以及工程實現（對前因后果和理論方面不感興趣的可以直接看這一小節的工程實現部分）

本章只是簡單介紹，如果想要深入了解在線凸優化（OCO），強烈推薦閱讀Elad Hazand的著作和Zinkevich的Paper

WM的理論證明可以參考Littlestone 94,Freund 99，雖然在上個世紀已經完成，但是將其理論拓展到一般的凸的函數還是在03年由Zinkevich完成的。

對於loss函數+正則化的結構風險最小化的優化問題（邏輯回歸也是這種形式）有兩種等價的描述形式，以1范數為例，分別是：

a、無約束優化形式的soft regularization formulation：

     

 

b、帶約束項的凸優化問題convex constraint formulation：

批量（batch）算法

 批量算法中每次迭代對全體訓練數據集進行計算（例如計算全局梯度），優點是精度和收斂還可以，缺點是無法有效處理大數據集（此時全局梯度計算代價太大），且沒法應用於數據流做在線學習。這里分無約束優化形式和約束優化（與上面問題描述可以對應起來）兩方面簡單介紹一下一些傳統批量算法。

• a、無約束優化形式：1、全局梯度下降，很常用的算法，就不細說了，每一步求一個目標函數的全局梯度，用非增學習率進行迭代；2、牛頓法（切線近似）、LBFGS（割線擬牛頓，用之前迭代結果近似Hessian黑塞矩陣的逆矩陣，BFGS似乎是幾個人名的首字母的簡稱）等方法。牛頓和擬牛頓等方法一般對於光滑的正則約束項（例如2范數）效果很好，據說是求解2范數約束的邏輯回歸類問題最好的方法，應用也比較廣，但是當目標函數帶L1非光滑、帶不可微點的約束項后，牛頓類方法比較無力，理論上需要做修改。感興趣的可以去查查無約束優化的相關數值計算的書，我也沒有更深入研究相關細節，這里不做重點關注。
• b、不等式約束凸優化形式：1、傳統的不等式約束優化算法內點法等；2、投影梯度下降（約束優化表示下），gt是subgradient，直觀含義是每步迭代后，迭代結果可能位於約束集合之外，然后取該迭代結果在約束凸集合上的投影作為新的迭代結果（第二個公式中那個符號標識向X的投影）：

在線算法

如上所述，批量算法有自身的局限性，而在線學習算法的特點是：每來一個訓練樣本，就用該樣本產生的loss和梯度對模型迭代一次，一個一個數據地進行訓練，因此可以處理大數據量訓練和在線訓練。常用的有在線梯度下降（OGD）和隨機梯度下降（SGD）等，本質思想是對上面【問題描述】中的未加和的單個數據的loss函數 L（w，zi)做梯度下降，因為每一步的方向並不是全局最優的，所以整體呈現出來的會是一個看似隨機的下降路線。典型迭代公式如下：

這里使用混合正則化項：，例如可能是1范數與2范數強凸項的混合（后面會看到其實很多都是這種混合正則化的格式，而且是有一定直觀含義的）。迭代公式中：gt是loss函數（單點的loss，未加和）的subgradient，與gt相加的那一項是混合正則化項中的第二項的梯度，投影集合C是約束空間（例如可能是1范數的約束空間），跟上面介紹的投影梯度下降類似的做法。

梯度下降類的方法的優點是精度確實不錯，但是不足相關paper主要提到兩點：

• 1、簡單的在線梯度下降很難產生真正稀疏的解，稀疏性在機器學習中是很看重的事情，尤其我們做工程應用，稀疏的特征會大大減少predict時的內存和復雜度。這一點其實很容易理解，說白了，即便加入L1范數（L1范數能引入稀疏解的簡單示例可以產看PRML那本書的第二章，我前面一篇blog的ppt里也大概提了），因為是浮點運算，訓練出的w向量也很難出現絕對的零。到這里，大家可能會想說，那還不容易，當計算出的w對應維度的值很小時，我們就強制置為零不就稀疏了么。對的，其實不少人就是這么做的，后面的Truncated Gradient和FOBOS都是類似思想的應用；
• 2、對於不可微點的迭代會存在一些問題，具體有什么問題，有一篇paper是這么說的：the iterates of the subgradient method are very rarely at the points of non-differentiability。我前后看了半天也沒看明白，有熟悉的同學可以指導一下。

三種做稀疏解的途徑

1）、簡單加入L1范數

• 局限如上面所提，a+b兩個float數很難絕對等於零，無法產生真正稀疏的特征權重

2）、在1范數的基礎上做截斷，最直觀沒技術含量的思路，那就設定一個閾值，做截斷來保證稀疏，可以結合L1范數

• 簡單截斷方法，每online訓練K個數據截斷一次，對OGD的迭代結果，每K步做一次截斷置零：

但是簡單截斷方法有問題：權重小，可能是確實是無用特征，還或者可能是該特征才剛被更新一次（例如訓練剛開始的階段、或者訓練數據中包含該特征的樣本數本來就很少），另外，簡單rounding技術太aggressive了，可能會破壞在線訓練算法的理論完備性。

• 簡單截斷基礎上，不太aggressive的Truncated gradient（TG） （09年的工作），其實后面的FOBOS也可以歸為這一類：

參數λ 和g_i控制模型稀疏，值越大越稀疏。 TG和簡單截斷法的區別在於 

3）、Black-box wrapper approaches：

• 黑盒的方法去除一些特征，然后重新訓練的看被消去的特征是否有效。
• 需要在數據集上對算法跑多次，所以不太實用

下面會提一下FOBOS（Forward-Backward Splitting method ，其實應該叫FOBAS的，歷史原因）以及RDA，因為后面的FTRL其實相當於綜合了這兩種算法的優點。

原始的OGD

優化函數：

其中x表示模型的參數，xi表示模型的第i個參數。

L(x)是邏輯回歸的似然函數的負對數。

最右邊那一項是L1正則項。

可以考慮加上L2正則來防止過擬合。

假設模型參數為x，第t個樣本的特征表示為vt

則樣本label為1的概率為：

該樣本對應的對數損失函數為：

令

    

則lt關於xt的求導過程如下：

上式即為第t個樣本的梯度g_t。

原始的OGD使用的迭代公式為：

其中et表示非增的學習率，也叫做步長，

OGD的算法比較簡單，不夠好，即使是加了L1正則化也可能不產生稀疏解。

FOBOS

FOBOS (Forward-Backward Splitting)是由John Duchi和Yoram Singer提出的^[11]。從全稱上來看，該方法應該叫FOBAS，但是由於一開始作者管這種方法叫FOLOS（Forward Looking Subgradients），為了減少讀者的困擾，作者干脆只修改一個字母，叫FOBOS。

該算法又稱為前后項算法，其權重的更新方式分為兩個部分 ,可以看作truncated gradient的一種特殊形式

基本思想：跟projected subgradient方法類似，不過將每一個數據的迭代過程，分解成一個經驗損失梯度下降迭代和一個最優化問題。分解出的第二個最優化問題，有兩項：第一項2范數那一項表示不能離第一步loss損失迭代結果太遠，第二項是正則化項，用來限定模型復雜度抑制過擬合和做稀疏化等。這個最優化問題有一些特殊的性質，從而保證了最終結果的稀疏性和理論上的完備，具體細節感興趣的可以查看對應paper。我這里更多關注直觀含義和工程實現，忽略理論方面的內容。

 

第一步為標准的梯度下降，然后將結果進行正則化的微調，滿足最速下降切保證了一定的稀疏性。其中ψ(W)ψ(W)是正則項可以對權重進行約束，常用的選擇為L1，則對應的方法為L1-FOBOS。

AOGD

此算法在網上極少相關文檔，還不能詳細介紹。

其迭代偽代碼為：

RDA

RDA（Regularized dual averaging），微軟10年的工作

• 非梯度下降類方法，屬於更加通用的一個primal-dual algorithmic schema的一個應用
• 克服了SGD類方法所欠缺的exploiting problem structure，especially for problems with explicit regularization。
• 能夠更好地在精度和稀疏性之間做trade-off

RDA（正則對偶平均）是微軟的研究成果，其權重更新策略為

權重更新包括三個部分第一部分線性加權，即對歷史梯度進行加權平均；第二部分正則項部分對特征進行稀疏化；第三部分嚴格遞增序列相當於額外的正則項。
實驗證明該方法能夠產生較好的稀疏和精度。

FTRL

主要用於CTR預測的在線訓練，成千上萬維度導致大量稀疏特征。一般希望模型參數更加稀疏，但是簡單的L1正則無法真正做到稀疏，一些梯度截斷方法（TG）的提出就是為了解決這個問題，在這其中FTRL是兼備精度和稀疏性的在線學習方法。FTRL的基本思想是將接近於0的梯度直接置零，計算時直接跳過以減少計算量。

在上文中我們從原理上定性比較了L1-FOBOS和L1-RDA在稀疏性上的表現。有實驗證明，L1-FOBOS這一類基於梯度下降的方法有比較高的精度，但是L1-RDA卻能在損失一定精度的情況下產生更好的稀疏性。那么這兩者的優點能不能在一個算法上體現出來？這就是FTRL要解決的問題。

FTRL（Follow the Regularized Leader）是由Google的H. Brendan McMahan在2010年提出的[1]，后來在2011年發表了一篇關於FTRL和AOGD、FOBOS、RDA比較的論文[2]，2013年又和Gary Holt, D. Sculley, Michael Young等人發表了一篇關於FTRL工程化實現的論文[3]。

本節我們首先從形式上將L1-FOBOS和L1-RDA進行統一，然后介紹從這種統一形式產生FTRL算法，以及介紹FTRL算法工程化實現的算法邏輯。

1.  L1-FOBOS和L1-RDA在形式上的統一

L1-FOBOS在形式上，每次迭代都可以表示為（這里我們令，是一個隨變化的非增正序列）：

 

 

把這兩個公式合並到一起，有：

 

 

通過這個公式很難直接求出的解析解，但是我們可以按維度將其分解為個獨立的最優化步驟：

由於與變量無關，因此上式可以等價於：

 

 

再將這N個獨立最優化子步驟合並，那么L1-FOBOS可以寫作：

 

 

而對於L1-RDA，我們可以寫作：

 

 

這里；如果令，，上面兩個式子可以寫作：

 

   公式(1)     公式(2)

 

需要注意，與論文[2]中的Table 1不同，我們並沒有將L1-FOBOS也寫成累加梯度的形式。

比較公式(1)和公式(2)，可以看出L1-FOBOS和L1-RDA的區別在於：(1) 前者對計算的是累加梯度以及L1正則項只考慮當前模的貢獻，而后者采用了累加的處理方式；(2) 前者的第三項限制的變化不能離已迭代過的解太遠，而后者則限制不能離點太遠。

2. FTRL算法原理

FTRL綜合考慮了FOBOS和RDA對於正則項和W限制的區別，其特征權重的更新公式為：

 

   公式(3)

 

注意，公式(3)中出現了L2正則項，在論文[2]的公式中並沒有這一項，但是在其2013年發表的FTRL工程化實現的論文[3]卻使用到了L2正則項。事實上該項的引入並不影響FRTL的稀疏性，后面的推導過程會顯示這一點。L2正則項的引入僅僅相當於對最優化過程多了一個約束，使得結果求解結果更加“平滑”。

公式(3)看上去很復雜，更新特征權重貌似非常困難的樣子。不妨將其進行改寫，將最后一項展開，等價於求下面這樣一個最優化問題：

 

 

上式中最后一項相對於來說是一個常數項，並且令，上式等價於：

 

 

針對特征權重的各個維度將其拆解成N個獨立的標量最小化問題：

 

 

到這里，我們遇到了與上一篇RDA中類似的優化問題，用相同的分析方法可以得到：

 

   公式(4)

 

從公式(4)可以看出，引入L2正則化並沒有對FTRL結果的稀疏性產生任何影響。

3.權重更新策略

工程優化：

1. 預測的內存方面：L1范式加策略，訓練結果w很稀疏，在用w做predict的時候節省了內存，很直觀

2. 訓練的內存方面：

1）在線丟棄訓練數據中很少出現的特征，即稀疏特征處理：

• Possion Inclusion(通過一定概率扔掉特征，對某一維度特征所來的訓練樣本，以p的概率接受並更新模型)
• BloomFilter Iclusion(通過碰撞優化特征，用bloom filter從概率上做某一特征出現k次才更新)

2）浮點數重新編碼

• 特征權重不需要用32bit或64bit的浮點數存儲，存儲浪費空間
• 16bit encoding，但是要注意處理rounding技術對regret帶來的影響

3）訓練若干相似model

• 對同一份數據序列，同時訓練多個相似的model
• 這些model有各自獨享的一些feature，也有一些共享的feature
• 出發點：有的特征維度可以是各個模型獨享的，而有的各個模型共享的特征，可以用同樣的數據訓練

4）Single Value Structure（據說有公司已經在實際中這么搞，大數據量下也能夠保證不錯的auc）

• 多個model公用一個feature存儲（例如放到cbase或redis中），各個model都更新這個共有的feature結構
• 對於某一個model，對於他所訓練的特征向量的某一維，直接計算一個迭代結果並與舊值做一個平均

5） 使用正負樣本的數目來計算梯度的和（所有的model具有同樣的N和P） 

６) Subsampling Training Data

• 在實際中，CTR遠小於50%，所以正樣本更加有價值。通過對訓練數據集進行subsampling，可以大大減小訓練數據集的大小
• 正樣本全部采（至少有一個廣告被點擊的query數據），負樣本使用一個比例r采樣（完全沒有廣告被點擊的query數據）。但是直接在這種采樣上進行訓練，會導致比較大的biased prediction
• 解決辦法：訓練的時候，對樣本再乘一個權重。權重直接乘到loss上面，從而梯度也會乘以這個權重，假設采樣概率為r，此時相當於每個樣本都有了權重w_t（正樣本權重為1，負樣本權重為1/r），同時相當於梯度也會乘上該系數，不會改變最終的損失。
• El(wt)=stwtl(wt)+(1−st)0=l(wt)，其中st=1/wt為該樣本被選中的概率
• 先采樣減少負樣本數目，在訓練的時候再用權重彌補負樣本，非常不錯的想法。

4. Per-Coordinate Learning Rates

前面介紹了FTRL的基本推導，但是這里還有一個問題是一直沒有被討論到的：關於學習率的選擇和計算。事實上在FTRL中，每個維度上的學習率都是單獨考慮的（Per-Coordinate Learning Rates）。

在一個標准的OGD里面使用的是一個全局的學習率策略，這個策略保證了學習率是一個正的非增長序列，對於每一個特征維度都是一樣的。

考慮特征維度的變化率：如果特征1比特征2的變化更快，那么在維度1上的學習率應該下降得更快。我們很容易就可以想到可以用某個維度上梯度分量來反映這種變化率。在FTRL中，維度i上的學習率是這樣計算的：

 

 公式(5)

 

由於，所以公式(4)中有，這里的和是需要輸入的參數，公式(4)中學習率寫成累加的形式，是為了方便理解后面FTRL的迭代計算邏輯。

5. FTRL算法邏輯與工程實現

大家對上面那一大坨前因后果和公式都不感興趣，ok，沒關系，google非常貼心地在13年給出了一篇工程性很強的paper，其實大部分公司使用FTRL的，根本不會關心上面那一大段東西，直接按着偽代碼寫，調調參，看結果很不錯就可以了。我們公司開始就是這么搞的，哈哈，不過人總是要有點兒好奇心的不是，深究一下前因后果和基本的理論公式感覺還是挺不同的。

邏輯回歸下的per-coordinate FTRL_Proximal的偽代碼如下，在公式表達的基礎上做了一些變換和實現上的trick，細節paper里有，大家在自己做實現的時候，可以在實際數據集上再並行加加速：

算法流程：

根據前面的學習，可以使用LR作為基本學習器，使用FTRL作為在線最優化的方法來獲取LR的權重系數，從而達到在不損失精度的前提下獲得稀疏解的目標。

工程化實現的幾個核心點是：

• 梯度計算代碼
• 損失函數計算代碼
• 權重更新計算代碼

算法代碼實現如下所示，這里只給出了核心部分代碼，主要做了以下優化:

• 修復原代碼更新梯度時候的邏輯錯誤。
• 支持pypy3加速，5.6G的訓練數據，使用Mac單機可以9分鍾跑完一個模型，一個epoch之后的logloss結果為0.3916；這對於刷競賽或者實際工程中模型訓練已經是比較理想的性能了。
• 使用16位浮點數保存權重數據，降低模型文件大小。
• 使用json保存模型非0參數w,z,n，進一步壓縮線上使用的時候占據的內存空間，可以進一步考慮壓縮模型文件大小，比如使用bson編碼+gzip壓縮后保存等。

from datetime import datetime
from csv import DictReader
from math import exp, log, sqrt
import gzip
import random
import json
import argparse


class FTRLProximal(object):
    """  FTRL Proximal engineer project with logistic regression  Reference:  https://static.googleusercontent.com/media/research.google.com/zh-CN//pubs/archive/41159.pdf   """

    def __init__(self, alpha, beta, L1, L2, D,
                 interaction=False, dropout=1.0,
                 dayfeature=True,
                 device_counters=False):

        # parameters
        self.alpha = alpha
        self.beta = beta
        self.L1 = L1
        self.L2 = L2
        self.dayfeature = dayfeature
        self.device_counters = device_counters

        # feature related parameters
        self.D = D
        self.interaction = interaction
        self.dropout = dropout

        # model
        self.n = [0.] * D
        self.z = [0.] * D
        self.w = [0.] * D

    def _indices(self, x):
        '''  A helper generator that yields the indices in x  The purpose of this generator is to make the following  code a bit cleaner when doing feature interaction.  '''

        for i in x:
            yield i

        if self.interaction:
            D = self.D
            L = len(x)
            for i in range(1, L):  # skip bias term, so we start at 1
                for j in range(i + 1, L):
                    # one-hot encode interactions with hash trick
                    yield abs(hash(str(x[i]) + '_' + str(x[j]))) % D

    def predict(self, x, dropped=None):
        """  use x and computed weight to get predict  :param x:  :param dropped:  :return:  """
        # wTx is the inner product of w and x
        wTx = 0.
        for j, i in enumerate(self._indices(x)):

            if dropped is not None and dropped[j]:
                continue

            wTx += self.w[i]

        if dropped is not None:
            wTx /= self.dropout

        # bounded sigmoid function, this is the probability estimation
        return 1. / (1. + exp(-max(min(wTx, 35.), -35.)))

    def update(self, x, y):
        """  update weight and coordinate learning rate based on x and y  :param x:  :param y:  :return:  """

        ind = [i for i in self._indices(x)]

        if self.dropout == 1:
            dropped = None
        else:
            dropped = [random.random() > self.dropout for i in range(0, len(ind))]

        p = self.predict(x, dropped)

        # gradient under logloss
        g = p - y

        # update z and n
        for j, i in enumerate(ind):

            # implement dropout as overfitting prevention
            if dropped is not None and dropped[j]:
                continue

            g_i = g * i
            sigma = (sqrt(self.n[i] + g_i * g_i) - sqrt(self.n[i])) / self.alpha
            self.z[i] += g_i - sigma * self.w[i]
            self.n[i] += g_i * g_i

            sign = -1. if self.z[i] < 0 else 1.  # get sign of z[i]

            # build w on the fly using z and n, hence the name - lazy weights -
            if sign * self.z[i] <= self.L1:
                # w[i] vanishes due to L1 regularization
                self.w[i] = 0.
            else:
                # apply prediction time L1, L2 regularization to z and get
                self.w[i] = (sign * self.L1 - self.z[i]) \
                            / ((self.beta + sqrt(self.n[i])) / self.alpha + self.L2)

    def save_model(self, save_file):
        """  保存weight數據到本地  :param save_file:  :return:  """
        with open(save_file, "w") as f:
            w = {k: v for k, v in enumerate(self.w) if v != 0}
            z = {k: v for k, v in enumerate(self.z) if v != 0}
            n = {k: v for k, v in enumerate(self.n) if v != 0}
            data = {
                'w': w,
                'z': z,
                'n': n
            }
            json.dump(data, f)

    def load_weight(self, model_file, D):
        """  loada weight data  :param model_file:  :return:  """
        with open(model_file, "r") as f:
            data = json.load(f)
            self.w = data.get('w', [0.] * D)
            self.z = data.get('z', [0.] * D)
            self.n = data.get('n', [0.] * D)

    @staticmethod
    def loss(y, y_pred):
        """  log loss for LR model  :param y:  :param y_pred:  :return:  """
        p = max(min(y_pred, 1. - 10e-15), 10e-15)
        return -log(p) if y == 1. else -log(1. - p)


def data(f_train, D, dayfilter=None, dayfeature=True, counters=False):
    ''' GENERATOR: Apply hash-trick to the original csv row  and for simplicity, we one-hot-encode everything   INPUT:  path: path to training or testing file  D: the max index that we can hash to   YIELDS:  ID: id of the instance, mainly useless  x: a list of hashed and one-hot-encoded 'indices'  we only need the index since all values are either 0 or 1  y: y = 1 if we have a click, else we have y = 0  '''

    device_ip_counter = {}
    device_id_counter = {}

    for t, row in enumerate(DictReader(f_train)):
        # process id
        ID = row['id']
        del row['id']

        # process clicks
        y = 0.
        if 'click' in row:
            if row['click'] == '1':
                y = 1.
            del row['click']

        # turn hour really into hour, it was originally YYMMDDHH

        date = row['hour'][0:6]
        row['hour'] = row['hour'][6:]

        if dayfilter != None and not date in dayfilter:
            continue

        if dayfeature:
            # extract date
            row['wd'] = str(int(date) % 7)
            row['wd_hour'] = "%s_%s" % (row['wd'], row['hour'])

        if counters:
            d_ip = row['device_ip']
            d_id = row["device_id"]
            try:
                device_ip_counter[d_ip] += 1
                device_id_counter[d_id] += 1
            except KeyError:
                device_ip_counter[d_ip] = 1
                device_id_counter[d_id] = 1
            row["ipc"] = str(min(device_ip_counter[d_ip], 8))
            row["idc"] = str(min(device_id_counter[d_id], 8))

        # build x
        x = [0]  # 0 is the index of the bias term
        for key in row:
            value = row[key]
            # one-hot encode everything with hash trick
            index = abs(hash(key + '_' + value)) % D
            x.append(index)
        yield t, ID, x, y

import tensorflow as tf
tf.compat.v1.train.FtrlOptimizer(learning_rate)

總結

實際的互聯網廣告應用需要的是快速地進行model的更新。為了保證快速的更新，訓練樣本是一條一條地過來的，每來一個樣本，model的參數對這個樣本進行一次迭代，從而保證了model的及時更新，這種方法叫做OGD（Online gradient descent）。

傳統Batch算法優點是精度和收斂還可以，缺點是無法有效處理大數據集（此時全局梯度計算代價太大），且沒法應用於數據流做在線學習

SGD存在的問題上面主要列了1）精度低；2）收斂慢；3）幾乎得不到稀疏解。其中對online learning最重要的問題是SGD很難得到需要的正則化設計的解，特別是幾乎得不到稀疏解

RDA，2010微軟提出，特點：相對FOBOS，在精度與稀疏性之間做平衡，其中實驗表明，在L1正則下，RDA比FOBOS可以更加有效地得到稀疏解。

 

與其他在線算法的迭代公式的對比（其實OGD如何一步步到類似形式的迭代公式的過程，限於時間，這里就不細說了，最后我會附一篇自己做分享會時做的ppt，里面有，感興趣的可以下載看看），不同的方法在這種統一的描述形式下，區別點僅在第二項和第三項的處理方式：

不同的方法按照統一的描述形式，區別點主要在第二項和第三項：

• 第一項：梯度或累積梯度；

• 第二項：L1正則化項的處理。

 

• 第三項：這個累積加和限定了新的迭代結果x不要離已迭代過的解太遠（也即FTRL-Proximal中proximal的含義），或者離0太遠（central），這一項其實也是low regret的需求

這小結做了在線最優化算法的整理和總結，沿着稀疏性的主線，先后介紹了簡單截斷法、TG、FOBOS、RDA以及FTRL。從類型上來看，簡單截斷法、TG、FOBOS屬於同一類，都是梯度下降類的算法,並且TG在特定條件可以轉換成簡單截斷法和FOBOS；RDA屬於簡單對偶平均的擴展應用；FTRL可以視作RDA和FOBOS的結合，同時具備二者的優點。目前來看，RDA和FTRL是最好的稀疏模型Online Training的算法。

對於Online模式的並行化計算，一方面可以參考ParallelSGD的思路，另一方面也可以借鑒batch模式下對高維向量點乘以及梯度分量並行計算的思路。總之，在理解算法原理的基礎上將計算步驟進行拆解，使得各節點能獨自無關地完成計算最后匯總結果即可。

三.相關注意問題

1.關於批量梯度下降的batch_size選擇問題

訓練集較小【<2000】：直接使用batch梯度下降，每次用全部的樣本進行梯度更新 
訓練集較大：batch_size一般設定為[64,512]之間，設置為2的n次方更符合電腦內存設置，代碼會運行快一些 
此外還要考慮GPU和CPU的存儲空間和訓練過程的波動問題 
 
batch_size越小，梯度的波動越大，正則化的效果也越強，自然訓練速度也會變慢，實驗時應該多選擇幾個batch_size進行實驗，以挑選出最優的模型。

2.關於批量梯度下降的weight_decay【似乎與L2正則有關系，待補充】

3.關於優化算法選擇的經驗之談

1. Adam在實際應用中效果良好，超過了其他的自適應技術。
2. 如果輸入數據集比較稀疏，SGD、NAG和動量項等方法可能效果不好。因此對於稀疏數據集，應該使用某種自適應學習率的方法，且另一好處為不需要人為調整學習率，使用默認參數就可能獲得最優值。
3. 如果想使訓練深層網絡模型快速收斂或所構建的神經網絡較為復雜，則應該使用Adam或其他自適應學習速率的方法，因為這些方法的實際效果更優。
4. SGD通常訓練時間更長，但是在好的初始化和學習率調度方案的情況下，結果更可靠。
5. Adadelta，RMSprop，Adam是比較相近的算法，在相似的情況下表現差不多。在想使用帶動量的RMSprop，或者Adam的地方，大多可以使用Nadam取得更好的效果。

4.訓練優化器的目的

加速收斂 2. 防止過擬合 3. 防止局部最優

5.選用優化器的目的

在構建神經網絡模型時，選擇出最佳的優化器，以便快速收斂並正確學習，同時調整內部參數，最大程度地最小化損失函數。

6.為什么神經網絡的訓練不采用二階優化方法 (如Newton, Quasi Newton)？

7. 優化SGD的其手段

 
 

import tensorflow as tf
import os
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
from matplotlib.animation import FuncAnimation
%matplotlib tk
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3'
plt.style.use('seaborn-white')

def func(x,y):
    return x**2-y**2

#制作畫圖用的網格數據
x1 = np.linspace(-1,1,100)
x2 = np.linspace(-1,2,100)
X1,X2 = np.meshgrid(x1,x2)
X = np.array(list(zip(X1.flatten(),X2.flatten())))
Y = np.array(list(map(lambda x:func(x[0],x[1]),X))).reshape(*X1.shape)
x = tf.Variable(initial_value=-0.5,dtype=tf.float32)
y = tf.Variable(initial_value=0.005,dtype=tf.float32)
z = tf.subtract(tf.square(x),tf.square(y))

def op_compare():
    sp = tf.Variable(1,dtype=tf.int64,trainable=False)
    sp_add = sp.assign(tf.add(sp,1))

    op_list = [
        tf.compat.v1.train.AdagradDAOptimizer(learning_rate=0.25,global_step=sp), #這個優化器一定要放在第一，否則會報錯，也搞不懂為什么。
        tf.compat.v1.train.GradientDescentOptimizer(learning_rate=0.25),
        tf.compat.v1.train.MomentumOptimizer(learning_rate=0.25,momentum=0.9),
        tf.compat.v1.train.MomentumOptimizer(learning_rate=0.25,momentum=0.9,use_nesterov=True,name='NAG'),
        tf.compat.v1.train.AdadeltaOptimizer(0.25),
        tf.compat.v1.train.RMSPropOptimizer(0.25),
        tf.compat.v1.train.AdamOptimizer(0.25)

    ]
    
    data_list = []
    with tf.compat.v1.Session() as sess:
        for i,init_op in enumerate(op_list):
            print('第{}個優化器'.format(i+1))
            op = init_op.minimize(z)
            sess.run(tf.compat.v1.global_variables_initializer())
            current_z = sess.run(z)
            diff = 2
            stop_e = 0.0001
            x_list = []
            print(sess.run(sp))
            res_x, res_y, res_z = sess.run([x,y,z])
            print(res_x,res_y)
            x_list.append([res_x,res_y,res_z])
            print(op.name)
            while diff>stop_e:
                sess.run([op,sp_add])
                res_x, res_y, res_z = sess.run([x,y,z])
                if res_z==np.inf or res_z==-np.inf or res_z<np.min(Y)-5:
                    break
                x_list.append([res_x,res_y,res_z])
                diff = abs(res_z-current_z)
                current_z = res_z
            data = np.array(x_list)
            data_list.append((op.name,data))
        return data_list

data_list = op_compare()

fig = plt.figure(figsize=(18,16))
ax = Axes3D(fig)
ax.plot_surface(X1,X2,Y,cmap=plt.cm.plasma)
count = 0
colors = ['r','g','b','y','purple','brown','pink']

line_list = [
    ax.plot([], [], [], color='{}'.format(c),marker='o', ls='-', markevery=[-1],markersize=5, animated=False,label='{}'.format(dt[0])) 
    for dt,c in zip(data_list,colors)
]   #重新設置曲線的值

def init():
    ax.set_ylim([-1, 2])
    ax.set_xlim([-1, 1])
    ax.set_zlim([np.min(Y),np.max(Y)])
    return line_list               #返回曲線

def update(n,data,lines):
#     print(n)
    for dt,ln in zip(data,lines):
        show_data = dt[1]
        id = n+1
        if n >= len(show_data)-1:
            id = -1
        xdata = show_data[:id,0]
        ydata = show_data[:id,1]
        zdata = show_data[:id,2]
        ln[0].set_data([xdata,ydata])
        ln[0].set_3d_properties(zdata)
    
    return lines
plt.legend()

#獲取數據的最大長度
record_len = []
for d in data_list:
    record_len.append(len(d[1]))
max_len = max(record_len)
ani = FuncAnimation(fig, update, frames=max_len, fargs=(data_list,line_list),
                    init_func=init,interval=100,)#這里的frames在調用update函數是會將frames作為實參傳遞給“n”

plt.show()

附錄