OWL-QN演算法

一、BFGS演算法

演算法思想如下：

Step1 取初始點 $x^{(0)}$ ，初始正定矩陣 $H_0$ ，允許誤差 $\epsilon>0$ ，令 $k=0$ ；

Step2 計算 $p^{(k)}=-H_k \nabla f(x^{(k)})$ ；

Step3 計算 $\alpha_k>0$ ，使得

$f(x^{(k)}+\alpha_kp^{(k)})={min}\limit_{\alpha \geq 0} f(x^{(k)}+\alpha p^{(k)})$ ；

Step4 令 $x^{(k+1)}=x^{(k)}+\alpha_k p^{(k)}$ ；

Step5 如果 $||\nabla f(x^{(k+1)})|| \leq \epsilon$ ，則取 $x^{(k+1)}$ 為近似最優解；否則轉下一步；

Step6 計算

$s_k=x^{(k+1)}-x^{(k)}$ ， $y_k=\nabla f(x^{(k+1)})-\nabla f(x^{(k)})$ ，

$H_{k+1}=H_k+\frac{1}{s_k^Ty_k}(1+\frac{y_k^TH_ky_k}{s_k^Ty_k})s_ks_k^T-\frac{1}{s_k^Ty_k}(s_ky_k^TH_k+H_ky_ks_k)$

令 $k=k+1$ ，轉Step2.

優點：

1、不用直接計算Hessian矩陣；

2、通過迭代的方式用一個近似矩陣代替Hessian矩陣的逆矩陣。

缺點：

1、矩陣儲存量為 $n^2$ ，因此維度很大時記憶體不可接受；

2、矩陣非稀疏會導致訓練速度慢。

二、L-BFGS演算法

針對BFGS的缺點，主要在於如何合理的估計出一個Hessian矩陣的逆矩陣，L-BFGS的基本思想是隻儲存最近的m次迭代資訊，從而大大降低資料儲存空間。對照BFGS，我重新整理一下用到的公式：

$\rho_k=\frac{1}{y_{k}^T s_k}$

$s_k=x_k-x_{k-1}$

$y_k=\nabla{f(x_k)}-\nabla{f(x_{k-1})}$

$V_k=I-\rho_{k}y_{k}s_{k}^T$

於是估計的Hessian矩陣逆矩陣如下：

$H_k=(I-\rho_{k-1}s_{k-1}y_{k-1}^T)H_{k-1}(I-\rho_{k-1}y_{k-1}s_{k-1}^T)+s_{k-1}\rho_{k-1}s_{k-1}^T$

$=V_{k-1}^TH_{k-1}V_{k-1}+ s_{k-1}\rho_{k-1}s_{k-1}^T$

把

$H_{k-1}=V_{k-2}^TH_{k-2}V_{k-2}+ s_{k-2}\rho_{k-2}s_{k-2}^T$

帶入上式，得：

$H_k=V_{k-1}^TV_{k-2}^TH_{k-2}V_{k-2}V_{k-1}+ V_{k-1}^Ts_{k-2}\rho_{k-2}s_{k-2}^T V_{k-1}+s_{k-1}\rho_{k-1}s_{k-1}^T$

假設當前迭代為k，只儲存最近的m次迭代資訊，（即：從k-m~k-1），依次帶入 $H$ ，得到：

公式1：

$H_k=(V_{k-1}^TV_{k-2}^T\ldots V_{k-m}^T) H_k^{0}(V_{k-m}\ldots V_{k-2}V_{k-1})$

$+ (V_{k-1}^TV_{k-2}^T\ldots V_{k-m+1}^T) S_{k-m}\rho_{k-m}S_{k-m}^T (V_{k-m}\ldots V_{k-2}V_{k-1})$

$+ (V_{k-1}^TV_{k-2}^T\ldots V_{k-m+2}^T) S_{k-m+1}\rho_{k-m+1}S_{k-m+1}^T (V_{k-m+1}\ldots V_{k-2}V_{k-1})$

$+\ldots$

$+V_{k-1}^T s_{k-2}\rho_{k-2} s_{k-2}^TV_{k-1}$

$+s_{k-1}\rho_{k-1}s_{k-1}^T$

演算法第二步表明了上面推導的最終目的：找到第k次迭代的可行方向，滿足：

$p_k=-H_k\nabla f(x_k)$

為了求可行方向p，有下面的:

two-loop recursion演算法

$q=\nabla f(x_k)$

$for (i=1 \ldots m) \quad \quad \quad do$

$\alpha_i=\rho_{k-i}s_{k-i}^Tq$

$q=q-\alpha_iy_{k-i}$

$end \quad \quad \quad for$

$r=H_k^{0}q$

$for (i=m \ldots 1) \quad \quad \quad do$

$\beta=\rho_{k-i}y_{k-i}^Tr$

$r=r+s_{k-i}(\alpha_i-\beta)$

$end \quad \quad \quad for$

$return \quad \quad \quad r$

該演算法的正確性推導：

1、令: $q_0=\nabla f(x_k)$ ，遞迴帶入q：

$q_i=q_{i-1}-\rho_{k-i}y_{k-i}s_{k-i}^Tq_{i-1}$

$=(I-\rho_{k-i}y_{k-i}s_{k-i}^T)q_{i-1}$

$=V_{k-i}q_{i-1}$

$=V_{k-i}V_{k-i+1}q_{i-2}$

$=\ldots$

$=V_{k-i}V_{k-i+1} \ldots V_{k-1} q_0$

$=V_{k-i}V_{k-i+1} \ldots V_{k-1} \nabla f(x_k)$

相應的：

$\alpha_i=\rho_{k-i}s_{k-i}^Tq_{i-1}$

$=\rho_{k-i}s_{k-i}^T V_{k-i+1}V_{k-i+2} \ldots V_{k-1} \nabla f(x_k)$

2、令： $r_{m+1}=H_{k-m}q=H_{k-m}V_{k-i}V_{k-i+1} \ldots V_{k-1} \nabla f(x_k)$

$r_i=r_{i+1}+s_{k-i}(\alpha_i-\beta) =r_{i+1}+s_{k-i}(\alpha_i-\rho_{k-i}y_{k-i}^Tr_{i+1})$

$=s_{k-i}\alpha_i+(I-s_{k-i}\rho_{k-i}y_{k-i}^T)r_{i+1}$

$=s_{k-i}\alpha_{i}+V_{k-i}^Tr_{i+1}$

於是:

$r_1=s_{k-1}\alpha_1+V_{k-1}^Tr_2 =s_{k-1}\rho_{k-1}s_{k-1}^T \nabla f(x_k)+V_{k-1}^Tr_2$

$=s_{k-1}\rho_{k-1}s_{k-1}^T \nabla f(x_k)+V_{k-1}^T(s_{k-2}\alpha_2+V_{k-2}^Tr_3)$

$=s_{k-1}\rho_{k-1}s_{k-1}^T \nabla f(x_k)+V_{k-1}^Ts_{k-2}\rho_{k-2}s_{k-2}^TV_{k-1}\nabla f(x_k)+V_{k-1}^T V_{k-2}^T r_3$

$=\ldots$

$=s_{k-1}\rho_{k-1}s_{k-1}^T \nabla f(x_k)$

$+V_{k-1}^T s_{k-2}\rho_{k-2} s_{k-2}^TV_{k-1} \nabla f(x_k)$

$+\ldots$

$+ (V_{k-1}^TV_{k-2}^T\ldots V_{k-m+2}^T) S_{k-m+1}\rho_{k-m+1}S_{k-m+1}^T (V_{k-m+1}\ldots V_{k-2}V_{k-1}) \nabla f(x_k)$

$+ (V_{k-1}^TV_{k-2}^T\ldots V_{k-m+1}^T) S_{k-m}\rho_{k-m}S_{k-m}^T (V_{k-m}\ldots V_{k-2}V_{k-1}) \nabla f(x_k)$

$+(V_{k-1}^TV_{k-2}^T\ldots V_{k-m}^T) H_{k-m}(V_{k-m}\ldots V_{k-2}V_{k-1}) \nabla f(x_k)$

這個two-loop recursion演算法的結果和公式1*初始梯度的形式完全一樣，這麼做的好處是：

1、只需要儲存 $s_{k-i}$ 、 $y_{k-i}$ （i=1~m）；

2、計算可行方向的時間複雜度從O(n*n)降低到了O(n*m)，當m遠小於n時為線性複雜度。

總結L-BFGS演算法的步驟如下：

Step1: 選初始點 $x_0$ ，允許誤差 $\epsilon >0$ ，儲存最近迭代次數m（一般取6）；

Step2: $k=0, \quad \quad \quad H_0=I ,\quad \quad \quad r=\nabla f(x_{0})$ ；

Step3: 如果 $||\nabla f(x_{k+1})||\leq \epsilon$ 則返回最優解 $x_{k+1}$ ，否則轉Step4；

Step4: 計算本次迭代的可行方向： $p_k=-r _k$ ；

Step5: 計算步長 $\alpha_k>0$ ，對下面式子進行一維搜尋：

$f(x_k+\alpha_kp_k)={min}\limits_{\alpha \geq 0} \quad \quad \quad f(x_k+\alpha p_k)$ ；

Step6: 更新權重x：

$x_{k+1}=x_k+\alpha_kp_k$ ；

Step7: if k > m

只保留最近m次的向量對，需要刪除( $s_{k-m},y_{k-m}$ )；

Step8: 計算並儲存：

$s_k=x_{k+1}-x_k$

$y_k=\nabla f(x_{k+1})-\nabla f(x_k)$ ；

Step9: 用two-loop recursion演算法求得：

$r_k=H_k\nabla f(x_k)$ ；

k=k+1，轉Step3。

需要注意的地方，每次迭代都需要一個 $H_{k-m}$ ，實踐當中被證明比較有效的取法為：

$H_k^0=\gamma_k I$

$\gamma_k=\frac{s_{k-1}^Ty_{k-1}}{{y_{k-1}^Ty_{k-1}}$

三、OWL-QN演算法

1、問題描述

對於類似於Logistic Regression這樣的Log-Linear模型，一般可以歸結為最小化下面這個問題：

$J(x)=l(x)+r(x)$

其中，第一項為loss function，用來衡量當訓練出現偏差時的損失，可以是任意可微凸函式（如果是非凸函式該演算法只保證找到區域性最優解），後者為regularization term，用來對模型空間進行限制，從而得到一個更“簡單”的模型。

根據對模型引數所服從的概率分佈的假設的不同，regularization term一般有：L2-norm（模型引數服從Gaussian分佈）；L1-norm（模型引數服從Laplace分佈）；以及其他分佈或組合形式。

L2-norm的形式類似於：

$J(x)=l(x)+C\sum\limit_i{x_i^2}$

L1-norm的形式類似於：

$J(x)=l(x)+C\sum\limit_i{|x_i|}$

L1-norm和L2-norm之間的一個最大區別在於前者可以產生稀疏解，這使它同時具有了特徵選擇的能力，此外，稀疏的特徵權重更具有解釋意義。

對於損失函式的選取就不在贅述，看兩幅圖：

圖1 - 紅色為Laplace Prior，黑色為Gaussian Prior

圖2 直觀解釋稀疏性的產生

對LR模型來說損失函式選取凸函式，那麼L2-norm的形式也是的凸函式，根據最優化理論，最優解滿足KKT條件，即有： $\nabla J(x^*)=0$ ，但是L1-norm的regularization term顯然不可微，怎麼辦呢？

2、Orthant-Wise Limited-memory Quasi-Newton

OWL-QN主要是針對L1-norm不可微提出的，它是基於這樣一個事實：任意給定一個維度象限，L1-norm 都是可微的，因為此時它是一個線性函式：

圖3 任意給定一個象限後的L1-norm

OWL-QN中使用了次梯度決定搜尋方向，凸函式不一定是光滑而處處可導的，但是它又符合類似梯度下降的性質，在多元函式中把這種梯度叫做次梯度，見維基百科http://en.wikipedia.org/wiki/Subderivative

舉個例子：

File:Subderivative illustration.png

圖4 次導數

對於定義域中的任何x₀，我們總可以作出一條直線，它通過點(x₀, f(x₀))，並且要麼接觸f的影象，要麼在它的下方。這條直線的斜率稱為函式的次導數，推廣到多元函式就叫做次梯度。

次導數及次微分：

凸函式f:I→R在點x₀的次導數，是實數c使得：

$f(x)-f(x_0)\ge c(x-x_0)$

對於所有I內的x。可以證明，在點x₀的次導數的集合是一個非空閉區間[a, b]，其中a和b是單側極限

$a=\lim_{x\to x_0^-}\frac{f(x)-f(x_0)}{x-x_0}$

$b=\lim_{x\to x_0^+}\frac{f(x)-f(x_0)}{x-x_0}$

它們一定存在，且滿足a ≤ b。所有次導數的集合[a, b]稱為函式f在x₀的次微分。

OWL-QN和傳統L-BFGS的不同之處在於：

1)、利用次梯度的概念推廣了梯度

定義了一個符合上述原則的虛梯度，求一維搜尋的可行方向時用虛梯度來代替L-BFGS中的梯度：

怎麼理解這個虛梯度呢？見下圖：

對於非光滑凸函式，那麼有這麼幾種情況：

圖5 $\partial_i^-f(x)>0$

圖6 $\partial_i^+f(x)<0$

圖7 otherwise

2)、一維搜尋要求不跨越象限

要求更新前權重與更新後權重同方向：

圖8 OWL-QN的一次迭代

總結OWL-QN的一次迭代過程：

–Find vector of steepest descent

–Choose sectant

–Find L-BFGS quadratic approximation

–Jump to minimum

–Project back onto sectant

–Update Hessian approximation using gradient of loss alone

最後OWL-QN演算法框架如下：

與L-BFGS相比，第一步用虛梯度代替梯度，第二、三步要求一維搜尋不跨象限，也就是迭代前的點與迭代後的點處於同一象限，第四步要求估計Hessian矩陣時依然使用loss function的梯度（因為L1-norm的存在與否不影響Hessian矩陣的估計）。

四、參考資料

1、Galen Andrew and Jianfeng Gao. 2007. 《Scalable training of L1-regularized log-linear models》. In Proceedings of ICML, pages 33–40.

2、http://freemind.pluskid.org/machine-learning/sparsity-and-some-basics-of-l1-regularization/#d20da8b6b2900b1772cb16581253a77032cec97e

3、http://research.microsoft.com/en-us/downloads/b1eb1016-1738-4bd5-83a9-370c9d498a03/default.aspx

from: http://www.cnblogs.com/vivounicorn/archive/2012/06/25/2561071.html