本文講解的是無約束優化中幾個常見的基於梯度的方法，主要有梯度下降與牛頓方法、BFGS 與 L-BFGS 算法。
梯度下降法是基於目標函數梯度的，算法的收斂速度是線性的，並且當問題是病態時或者問題規模較大時，收斂速度尤其慢（幾乎不適用）；
牛頓法是基於目標函數的二階導數（Hesse 矩陣）的，其收斂速度較快，迭代次數較少，尤其是在最優值附近時，收斂速度是二次的。但牛頓法的問題在於當海森矩陣稠密時，每次迭代的計算量比較大，因為每次都會計算目標函數的海森矩陣的逆，這樣一來，當數據維度較高時，不僅計算量大（有時大到不可計算），而且需要的存儲空間也多，因此牛頓法在面對海量數據時由於每一步迭代的開銷巨大而變得不適用。
擬牛頓法是在牛頓法的基礎上引入了海森矩陣的近似矩陣，避免每次迭代都要計算海森矩陣的逆，擬牛頓法的收斂速度介於梯度下降法和牛頓法之間，是超線性的。擬牛頓法的問題也是當問題規模很大時，近似矩陣變得很稠密，在計算和存儲上也有很大的開銷，因此變得不實用。
另外需要注意的是，牛頓法在每次迭代時不能總是保證海森矩陣是正定的，一旦海森矩陣不是正定的，優化方向就會“跑偏”，從而使得牛頓法失效，也說明了牛頓法的魯棒性較差。擬牛頓法用海森矩陣的逆矩陣來替代海森矩陣，雖然每次迭代不能保證是最優的優化方向，但是近似矩陣始終是正定的，因此算法總是朝着最優值的方向在搜索。這些方法均是基於迭代的方法，通過找到序列 $x_1,x_2,….x_k,…$ 來一步步極小化目標函數，達到求取極小值的目的。
大概了解了各種算法的特點之后，接下來步入正文，對於無約束優化問題，變量 $x \in \mathbb{R}^n$ ,目標函數為：
\[\min_xf(x)\]
泰勒級數
基於梯度的方法都會涉及泰勒級數問題，這里簡單介紹一下，泰勒級數就是說函數 $f(x)$ 在點 $x_0$ 的鄰域內具有 $n+1$ 階導數，則該鄰域內 $f(x)$ 可展開為 $n$ 階泰勒級數為：
\[f(x) = f(x_0) + \nabla f(x_0)(x-x_0)+\frac{\nabla^2f(x_0)}{2!}(x-x_0)^2+…+\frac{\nabla^nf(x_0)}{n!}(x-x_0)^n\]
梯度下降法
梯度下降法是一種通過迭代來求取極值的方法，它使用目標函數的一階導數信息，每次迭代取目標函數值下降最快的方向作為搜索方向，目標函數 $f(x)$ 在什么方向下降最快呢？這里給出一個證明，假設當前為第 $k$ 次迭,在自變量取值 $x_k$ 處進行泰勒展開可得：
\[f(x) = f(x_k) + \nabla f(x_k)(x-x_k)\]
這里嚴格意義上應該取 “ $\approx$ ” 因為只取了一階泰勒，之后的牛頓法，展開為二階泰勒也是近似相等的關系。現在 $f(x)$ 即為下一步的取值的函數，將 $x – x_k $ 即為下一步要迭代的方向，記做 $ag$，意思是說 $x$ 沿單位向量 $g \in \mathbb{R}^n$ 的方向迭代，$\alpha$ 為步長參數記 ，代表走多大的一步，現在一階泰勒可寫作：
\[f(x_k) – f(x) = – \alpha \nabla f(x_k) \cdot g \]
每次迭代是為了使目標函數盡可能減小，即使得 $f(x_k) – f(x)$ 盡可能大，忽略步長參數 $\alpha$， 可得極大化 $f(x_k) – f(x)$ 等價於：
\[ \min \nabla f(x_k) \cdot g\]
$\nabla f(x) \cdot g$ 代表兩向量的乘積，很明顯，兩向量方向相反時便能取得最小值，所以要想使得迭代后以目標函數以最快速度下降，迭代的方向 $g$ 應該是梯度的反方向， $x$ 每次沿着 $-\nabla f(x)$ 的方向前進便能得到極小值，綜上梯度下降的算法如下：
\[x^{k+1} := x^{k} – \alpha \nabla f(x_k)\]
負梯度是下降最快的方向，梯度是上升最快的方向，證明就是上述的很簡單的一個泰勒展開。
牛頓方法
牛頓方法也是一種迭代的方法，不同於梯度下降，該方法引入了二階導數信息，假設當前迭代到第 $k$ 次，將目標函數在自變量 $x_k$ 處展開為二階泰勒級數：
\[ f(x) = f(x_k) + \nabla f(x_k) (x-x_k)+ \frac{1}{2} (x –x_k)^T\nabla^2 f(x_k )(x-x_k)\]
在這個關於 $x$ 的二次函數里，另其導數得 $0$ ，得到的點即可作為下一次自變量的值 $x_{k+1}$，人們常說的牛頓方法是曲面擬合，這個曲面正是關於 $x$ 的二次函數構成的，而在梯度法中由於一階泰勒只用到 $x$ 的一次函數，所以梯度法是平面擬合，下圖展示了牛頓法與梯度法的區別：
接下來對 $f(x)$ 兩端同時對 $x$ 求導，另導數得零求得的極值點就是下一次迭代的取值 $x_{k+1}$：
\[\nabla f(x) = \nabla f(x_k) + \nabla^2 f(x_k )(x-x_k) = 0 \]
公式里的 $\nabla f$ 為梯度 ,$\nabla^2 f$ 為 Hesse 矩陣，方便起見，將 $\nabla f$ 記做 $g$ ,$\nabla^2 f$ 記做 $H$ ， $g$ 與 $H$ 的形式分別如下：
\[ g = \nabla f = \begin{bmatrix} \frac{\partial f }{\partial x_1} \\ \frac{\partial f }{\partial x_2} \\ \vdots\\\frac{\partial f }{\partial x_n} \\ \end{bmatrix} \ \ \ \
H = \nabla^2 f = \begin{bmatrix}
\frac{\partial^2 f }{\partial^2 x_1}& \frac{\partial^2 f }{\partial x_1 x_2} & \cdots &\frac{\partial ^2 f }{\partial x_1x_n} \\
\frac{\partial^2 f }{\partial^2 x_1x_2}& \frac{\partial^2 f }{\partial x_2^2} &\cdots &\frac{\partial ^2 f }{\partial x_2x_n} \\
\vdots & \vdots & \ddots & \vdots \\
\frac{\partial^2 f }{\partial^2 x_nx_1}& \frac{\partial^2 f }{\partial x_n x_2} &\cdots &\frac{\partial ^2 f }{\partial x_n^2} \\
\end{bmatrix} \]
采用符號 $g$ 與 $H$ 表示：
\[g_k + H_k(x-x_k) = 0\]
綜上於是 $x$ 在下一次迭代的取值 ：
\[x_{k+1} = x_k -H^{-1}_k g_k \]
這里 $-H^{-1}_k g_k$ 即為搜索方向向量，函數值需要在此方向下降，所以需要與梯度 $g_k$ 反向 ，即 $-g_kH^{-1}_k g_k<0$ ，顯而易見需要 Hesse 矩陣是正定的。如果這一條件不滿足，可能導致目標函數不一定下降，從而牛頓法不收斂。比如說初始點 $x_0$ 離極小值點比較遠，就可能造成其 Hesse 矩陣不是正定的，甚至是正定，但目標函數值沒有下降。對於二次型函數，牛頓法一次迭代便可達到極值點，但對於非二次型，牛頓法缺少步長因子，所以迭代可能導致發散或者不穩定的現象，所以在迭代時，需要在每一步迭代時給出一個步長參數 $\lambda_k$：
\[ \lambda_k = arg \min_{\lambda} f(x_k – \lambda H_k^{-1}g_k)\]
綜上，總結一下牛頓法的特點：
- 牛頓法每次計算都需要計算目標函數的 Hesse 矩陣，計算量非常大，且要求 Hesse矩陣正定；
- 梯度下降法是一階收斂，牛頓法是二階收斂，由於二階導 $\Delta (x_k)$ 可能為0甚至為負（非正定），則迭代后的 $x_{k+1}$ 可能不是下降方向（如果為0，則無法計算）。
- 牛頓法是用一個二次曲面去擬合當前所處位置的局部曲面，而梯度下降法是用一個平面去擬合當前的局部平面，所以牛頓法選擇的下降路經會更符合真實的最優下降路經；
下圖來自 wiki ，紅色的為牛頓法，搜索方向明顯優於地圖下降法綠線：
擬牛頓法
擬牛頓法可以克服牛頓法計算量大的缺點，不在計算目標函數的 Hesse 矩陣，而是構造一個近似 Hesse 矩陣的對稱正定矩陣，根據近似矩陣來優化目標函數，不同的近似構造 Hesse 的方法決定了不同的擬牛頓法，構造 Hesse 矩陣是需要滿足擬牛頓條件的，擬牛頓條件是這樣求得的，首先將 $f(x)$ 在 $x_{k+1}$ 處做二階泰勒展開：
\[ f(x) = f(x_{k+1}) + \nabla f(x_{k+1}) (x-x_{k+1})+ \frac{1}{2} (x –x_{k+1})^T\nabla^2 f(x_{k+1} )(x-x_{k+1})\]
兩邊同時對 $x_{k+1}$求導可得：
\[g = g_{k+1}+ H_{k+1}(x -x_{k+1}) \]
令 $x = x_k$ ，整理可得：
\[ g_{k+1} – g_k = H_{k+1} (x_{k+1} – x_k)\]
這個便是擬牛頓條件了，迭代過程中對 $H_{k+1}$ 做出約束，根據約束構造一個近似矩陣 $B_{k+1}$ ,來模擬 Hesse 矩陣就可以了，為了簡便起見，引入記號 $s_k$ 與 $y_k$ ：
\[ s_k = x_{k+1} –x_k , y_k = g_{k+1} –g _k \]
因此可得擬牛頓條件為：
\[ y_k = B_{k+1} \cdot s_k \]
因為牛頓法中的迭代方向為 $-H^{-1} \cdot g$，所以另 D_{k+1} = H_{k+1}^{-1}，擬牛頓條件還可以寫作：
\[s_k = D_{k+1} \cdot y_k \]
BFGS
BFGS 是一種擬牛頓方法，通過迭代構建近似 Hesse 矩陣，省去了求解 Hesse 的復雜的步驟，而且 BFGS 構造出來的近似 Hesse 矩陣一定是正定的，這完全克服了牛頓法的缺陷，雖然搜索方向不一定最優，但始終朝着最優的方向前進的。首先初始化 Hesse 矩陣 $B_0 = I $ ，接下來每次迭代對矩陣 $B_k$ 進行更新即可：
\[ B_{k+1} = B_k+ \Delta B_k , \ k = 1,2,…\]
迭代構建近似矩陣的關鍵是矩陣 $\Delta B_k$ 的構造了，將其寫作：
\[\Delta B_k = \alpha uu^T + \beta vv^T\]
這里的向量 u 和 v 是待定的，知道了這兩個向量，就可以構造構造 $\Delta B_k$ 了，且這樣構造出的矩陣是對稱的，根據擬牛頓條件：
\begin{aligned}
y_k &= B_{k+1} s_k \\
&= (B_k + \Delta B_k)s_k \\
&= (B_k + \alpha uu^T + \beta vv^T)s_k \\
&= B_k s_k + (\alpha u^Ts_k) \cdot u+ (\beta v^Ts_k) \cdot v
\end{aligned}
這里 $\alpha u^Ts_k$ 與 $\beta v^Ts_k$ 均為實數，代表了 在 $u$ 與 $v$ 方向的拉伸程度，為了計算簡單，做如下賦值運算：
\[\alpha u^Ts_k = 1 , \ \beta v^Ts_k = –1\]
代入上式便可得：
\[u - v = y_k – B_ks_k\]
這就得到得到了 $u$ 與 $v$ 的一個近似:
\[ u = y_k , \ v = B_k s_k\]
繼而求 $\alpha$ 與 $\beta$ 的值，下式推倒中用了$B_k$ 對稱的特征：
\[\alpha = \frac{1}{y^Ts_k}, \beta= -\frac{1}{(B_ks_k)^Ts_k} = -\frac{1}{s_k^TB_ks_k}\]
$\alpha$ 、 $\beta$ 、 $u$ 與 $v$都求得后，便得到了 $\Delta B_k$ 的更新公式：
\[\Delta B_k = \frac{y_ky_k^T}{y_k^Ts_k} – \frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}\]
因此迭代計算的公司為：
\[B_{k+1} = B_k +\frac{y_ky_k^T}{y_k^Ts_k} – \frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}\]
由與牛頓法的方向是 $–H^{-1}_kg_k$ 的，所以最好可以直接計算出 $B_k^{-1}$ ，這樣就不用再進行求逆運算了，直接根據Sherman-Morrison 公式：可得關於矩陣B 的逆的更新方式：
\[B^{-1}_{k+1} = B^{-1}_k + \left (\frac{1}{s_k^Ty_k}+\frac{y_k^TB_k^{-1}y_k}{(s_k^Ty_k)^2} \right )s_ks_k^T - \frac{1}{s_k^Ty_k} \left (B_k^{-1}y_ks_k^T + s_ky_k^TB^{-1}_k \right)\]
$B^{-1}_k$ 不就是 $D_k$ 嗎，這里用 $D_k$ 來表示，給出最終的 BFGS 算法：
1. 給出 $x_0$ 與迭代停止條件 ，另 $D_0 = I, \ k = 0$;
2. $for$ $k = 1,2,…K$, $do$:
2.1 牛頓法的搜索方向為 $d_k = – D_kg_k$
2.2 依照牛頓方法計算步長因子 $\lambda_k$:\[s_k = -\lambda_k d_k, \ \ x_{k+1} = x_k + s_k\]
2.3 滿足條件則停止迭代
2.4 計算 $y_k = g_{k+1} – g_k$
2.6 更新擬Hesse 矩陣：
\[D_{k+1} = D_k + \left (\frac{1}{s_k^Ty_k}+\frac{y_k^T D_ky_k}{(s_k^Ty_k)^2} \right )s_ks_k^T - \frac{1}{s_k^Ty_k}(D_ky_ks_k^T + s_ky_k^T D_k )\]
這就是全部的 BFGS 算法了
停止條件為人為設定，可設定為兩次迭代目標函數差的閾值或者梯度差的閾值，或者梯度本身作為閾值。
L-BFGS
工業中實用的擬牛頓法的便是 L-BFGS 了，以前一直只聽說過它的大名，其實就是一種擬牛頓方法而已，這里也不細看了，對於近似 Hesse 矩陣 $D_k$ ：
\[D_{k+1} = D_k + \left (\frac{1}{s_k^Ty_k}+\frac{y_k^T D_ky_k}{(s_k^Ty_k)^2} \right )s_ks_k^T - \frac{1}{s_k^Ty_k}(D_ky_ks_k^T + s_ky_k^T D_k )\]
而是存儲向量序 $s_k,y_k$,而且向量序列也不是都存，而是存最近的 $m$ 次的， $m$ 為人工指定，計算 $D_k$ 時，只用最新的 $m$ 個向量模擬計算即可。在第 k 次迭代，算法求得了 $x_k$ ，並且保存的曲率信息為 ${s_i, y_i}_{k-m}^{k-1}$。為了得到 $H_k$，算法每次迭代均需選擇一個初始的矩陣 $H_K^0$，這是不同於 BFGS 算法的一個地方，接下來只用最近的 m 個向量對該初始矩陣進行修正，實踐中 $H_k^0$ 的設定通常如下：
\begin{aligned}
H_k^0 &=r_kI \\
r_k &=\frac{s{k-1}^Ty_{k-1}}{y_{k-1}^Ty_{k-1}}
\end{aligned}
其中 $r_k$ 表示比例系數，它利用最近一次的曲率信息來估計真實海森矩陣的大小，這就使得當前步的搜索方向較為理想，而不至於跑得“太偏”，這樣就省去了步長搜索的步驟，節省了時間。在L-BFGS算法中，通過保存最近 m 次的曲率信息來更新近似矩陣的這種方法在實踐中是很有效的，雖然 L-BFGS 算法是線性收斂，但是每次迭代的開銷非常小，因此 L-BFGS 算法執行速度還是很快的，而且由於每一步迭代都能保證近似矩陣的正定，因此算法的魯棒性還是很強的。
總結下 BFGS 與 L-BFGS 的： BFGS算法在運行的時候，每一步迭代都需要保存一個 $n \times n$ 的矩陣，現在很多機器學習問題都是高維的，當 $n$ 很大的時候，這個矩陣占用的內存是非常驚人的，並且所需的計算量也是很大的，這使得傳統的 BFGS 算法變得非常不適用。而 L-BFGS 則是很對這個問題的改進版，從上面所說可知，BFGS 算法是通過曲率信息 $(sk, yk)$ 來修正 $H_k$ 從而得到 $Hk+1$ ，L-BFGS 算法的主要思路是：算法僅僅保存最近 $m$ 次迭代的曲率信息來計算 $H_{k+1}$ 。這樣，我們所需的存儲空間就從 $n \times n$ 變成了 $2 m \times n$ 而通常情況下 $m<<n$ 。
參考文獻：
http://tangxman.github.io/2015/11/19/optimize-newton/
http://www.cnblogs.com/richqian/p/4535550.html
http://www.cnblogs.com/zhangchaoyang/articles/2600491.html
http://mlworks.cn/posts/introduction-to-l-bfgs/
http://blog.csdn.net/itplus/article/details/21896981
http://blog.csdn.net/henryczj/article/details/41542049
http://www.voidcn.com/blog/dadouyawp/article/p-4204403.html 很棒的文章