梯度下降法原理與模擬分析||系列（1）

1 引言

梯度下降法（Gradient Descent）也稱為最速下降法（Steepest Descent），是法國數學家奧古斯丁·路易·柯西 (Augustin Louis Cauchy) 於1847年提出來，它是最優化方法中最經典和最簡單的一階方法之一。梯度下降法由於其較低的複雜度和簡單的操作而在很多領域得到廣泛研究和應用，如機器學習。由梯度下降法衍生了許多其他演算法，如次梯度下降法，近端梯度下降法，隨機梯度下降法，回溯梯度發，動量加速梯度法等等。本文只介紹最基礎的梯度下降法原理和理論分析，與此同時，通過模擬來說明梯度下降法的優勢和缺陷。其他重要的梯度下降衍生方法會持續更新，敬請關注。

2 梯度下降法原理

2.1 偏導數，方向導數和梯度

在直角座標系下，標量函式 \(f:\mathbb{R}^{n}\mapsto \mathbb{R}\) 的梯度 \(\nabla f\) 定義為：

\[\nabla f = \frac{{\partial f}}{{\partial x}}{\bf{i}} + \frac{{\partial f}}{{\partial y}}{\bf{j}} + \frac{{\partial f}}{{\partial z}}{\bf{k}} \]

其中 \(\bf{i}\)，\(\bf{j}\)，\(\bf{k}\) 表示在三個維度的單位方向向量。\(\frac{{\partial f}}{{\partial x}}\)，\(\frac{{\partial f}}{{\partial y}}\)，\(\frac{{\partial f}}{{\partial z}}\) 為對應的偏導數。

• 方向導數：函式在某點處有無數條切線，沿著這些方向的斜率組成了方向導數，每一條切線都代表一個變化的方向。
• 偏導數：函式在一個點處有無數個導數，而延著座標軸方向對應的導數叫偏導數。
• 梯度：是一個向量，它是的方向是最大方向導數所對應的方向。

梯度下降法就是以梯度為搜尋方向的迭代優化演算法

圖2. 梯度下降法迭代過程

2.2 梯度下降法描述

對於無約束優化問題：

\[\mathop {\arg \min }\limits_{{\bf{x}} \in {\mathbb{R}^n}} f({\bf{x}}) \]

函式 \(f\) 可以是凸的，也可以是非凸的，在高維非凸的情況下，解這樣一個優化問題可能很難計算解析解，而我們可以使用數值迭代法來求得一個區域性最優解。這樣的迭代方法有很多，比如本文講的梯度下降法，牛頓迭代法，共軛梯度法等，其中梯度下降法利用了梯度資訊，屬於一階方法，牛頓法利用了海森矩陣，屬於二階方法，而共軛梯度法是介於這二者之間的方法。

圖1. 複雜函式影像

我們先來看看梯度下降法是什麼樣子的，梯度下降法的迭代公式為：

\[{{\bf{x}}^{k + 1}} = {{\bf{x}}^k} + {\alpha _k}{{\bf{d}}^k} \]

式 (2) 中 \({{\bf{x}}^k}\) 表示第 \(k\) 次迭代時的位置，\(\alpha_k\) 表示第 \(k\) 次迭代的步長，\(\bf{d}_k\) 表示第 \(k\) 次迭代的尋優方向。梯度下降法就是在給定初始點 \(\bf{x}_0\) 後，通過不斷沿著尋優方向迭代找到區域性最優值的過程。

那麼梯度下降法中的步長和方向怎麼確定呢？
一般來說，步長的選擇有多種方式，可以是固定的，也可以是隨著迭代過程而變化的
尋優方向則是 \(f(\bf{x})\) 在點 \(\bf{x}_k\) 處的梯度反方向，

2.3 梯度下降法的矩陣形式

以求解線性方程 \({\bf{y}} = {\bf{Ax}}\) 為例，
假設目標函式為

\[f({\bf{x}}) = \frac{1}{2}\left\| {{\bf{Ax}} - {\bf{y}}} \right\|_2^2 \]

我們知道目標函式的梯度就是梯度下降法的尋優方向，所以我們可以先求目標函式的導數：

\[{\bf{g}} = {{\bf{A}}^T}({\bf{Ax}} - {\bf{y}}) \]

在 \({\bf{x}}_{k}\) 處，將 \({{\bf{x}}_{k + 1}} = {{\bf{x}}_k} + {\alpha _k}{{\bf{g}}_k}\) 代入上式中，得到：

\[f({{\bf{x}}_{k + 1}}) = \frac{1}{2}\left\| {{\bf{A}}({{\bf{x}}_k} + {\alpha _k}{{\bf{g}}_k}) - {\bf{y}}} \right\|_2^2 \]

接下來，為了求 \({\bf{x}}_{k}\) 處的步長，我們介紹一種精確線搜尋（Exact Line Search）方法，
精確線搜尋方法的原理是沿著 \({{\bf{x}}_k} + {\alpha _k}{{\bf{g}}_k}\) 方向最小化目標函式，從而得到此時的 \({\alpha _k}\)。也就是求解這樣一個問題：

\[{\alpha _k} = \mathop {\arg \min }\limits_\alpha f\left( {{{\bf{x}}_k} + {\alpha _k}{{\bf{g}}_k}} \right) \]

求解這個優化問題，我們只需對 \({\alpha _k}\) 求導，並且令結果等於0，也就是：

\[\frac{{\partial f({{\bf{x}}_{k + 1}})}}{{\partial {\alpha _k}}} = 0 \]

\[{\left( {{\bf{A}}{{\bf{g}}_k}} \right)^T}\left( {{\bf{A}}({{\bf{x}}_k} + {\alpha _k}{{\bf{g}}_k}) - {\bf{y}}} \right) = 0 \]

\[{\alpha _k}{\left( {{\bf{A}}{{\bf{g}}_k}} \right)^T}{\bf{A}}{{\bf{g}}_k} + {\bf{g}}_k^T{{\bf{A}}^T}\left( {{\bf{A}}{{\bf{x}}_k} - {\bf{y}}} \right) = 0 \]

\[{\alpha _k}{\left( {{\bf{A}}{{\bf{g}}_k}} \right)^T}{\bf{A}}{{\bf{g}}_k} = {\bf{g}}_k^T{{\bf{g}}_k} \]

\[{\alpha _k} = \frac{{{\bf{g}}_k^T{{\bf{g}}_k}}}{{{{\left( {{\bf{A}}{{\bf{g}}_k}} \right)}^T}{\bf{A}}{{\bf{g}}_k}}} \]

當然，還有其他的一些步長選擇方法，比如回溯線搜尋方法等，在這裡就不再展開。
尋優方向和步長都得到了，那麼矩陣形式的梯度下降法的迭代公式可以總結為：

\[{{\bf{x}}_{k + 1}} = {{\bf{x}}_k} + {\alpha _k}{{\bf{A}}^T}({\bf{A}}{{\bf{x}}_k} - {\bf{y}}) \]

可以證明，採用精確線性搜尋的最速下降法的收斂速度是線性的。對於極小化正定二次函式，梯度下降法產生的序列滿足

\[\frac{{f({x_{k + 1}}) - f({x^*})}}{{f({x_k}) - f({x^*})}} \le {\left( {\frac{{{\lambda _1} - {\lambda _n}}}{{{\lambda _1} + {\lambda _n}}}} \right)^2} = {\left( {\frac{{\kappa - 1}}{{\kappa + 1}}} \right)^2} \]

這裡 \({{\lambda _1}}\) 和 \({{\lambda _1}}\) 分別表示矩陣 \({{\bf{A}}^T}{\bf{A}}\) 的最大和最小特徵值，\(\kappa = \frac{{{\lambda _1}}}{{{\lambda _n}}}\) 表示矩陣 \({{\bf{A}}^T}{\bf{A}}\) 的條件數。可以看出梯度下降法的收斂速度是與對稱矩陣 \({{\bf{A}}^T}{\bf{A}}\) 的條件數有關的，條件數越小收斂速度越快。如果矩陣的條件數遠大於1，則說明矩陣是ill-conditioned，反之，如果矩陣的條件數等於1，則說明矩陣是well-conditioned。

3 收斂性分析

如果函式 \(f\) 在定義域內是利普希茲連續（Lipschitz continuous）的，則有：

\[\left\| {f(x) - f(y)} \right\| \le L\left\| {x - y} \right\| \]

如果函式 \(f\) 的導數符合Lipschitz continuous，那麼有：

\[\left\| {\nabla f(x) - \nabla f(y)} \right\| \le L\left\| {x - y} \right\| \]

其中 \(L>0\) 稱為利普希茲常數。

在梯度滿足利普希茲連續條件下，梯度下降法的收斂性如下：

步長為 \(\alpha \le \frac{1}{L}\) 的梯度下降法滿足：

\[f({x_k}) - f({x^*}) \le \frac{{{{\left\| {{x_0} - {x^*}} \right\|}^2}}}{{2\alpha k}} \]

也就是說梯度下降法的收斂速度為 \(O\left( {\frac{1}{k}} \right)\)，只需要 \(O\left( {\frac{1}{\varepsilon }} \right)\) 次迭代即可達到 $f({x_k}) - f({x^*}) \le \varepsilon $。

上述結論的證明見參考文獻[1]，這裡給一個簡單的描述。
如果 \(\nabla f(x)\) 滿足利普希茲連續條件，那麼下式成立

\[f(y) \le f(x) + \nabla f{(x)^T}(y - x) + \frac{L}{2}{\left\| {y - x} \right\|^2} \]

將 \(y = x - \alpha \nabla f(x)\) 代入可得

\[f(y) \le f(x) - (1 - \frac{{L\alpha }}{2})\alpha {\left\| {\nabla f(x)} \right\|^2} \]

令 \(\tilde x = x - \alpha \nabla f(x)\)，且 \(0 < \alpha \le \frac{1}{L}\)，有

\[\begin{array}{l} f(\tilde x) \le f({x^*}) + \nabla f{(x)^T}(y - {x^*}) - \frac{\alpha }{2}{\left\| {\nabla f(x)} \right\|^2} \\ = f({x^*}) + \frac{1}{{2\alpha }}({\left\| {x - {x^*}} \right\|^2} - {\left\| {\tilde x - {x^*}} \right\|^2}) \\ \end{array}\]

將所有迭代求和

\[\begin{array}{l} \sum\limits_{i = 1}^k {\left( {f({x_i}) - f({x^*})} \right)} \le \frac{1}{{2\alpha }}\left( {{{\left\| {{x_0} - {x^*}} \right\|}^2} - {{\left\| {{x_k} - {x^*}} \right\|}^2}} \right) \\ \le \frac{1}{{2\alpha }}{\left\| {{x_0} - {x^*}} \right\|^2} \\ \end{array}\]

因為 \(f({x_k})\) 是非增的，所以有

\[f({x_i}) - f({x^*}) \le \frac{1}{k}\sum\limits_{i = 1}^k {\left( {f({x_i}) - f({x^*})} \right)} \le \frac{{{{\left\| {{x_0} - {x^*}} \right\|}^2}}}{{2\alpha k}} \]

4 梯度下降法模擬

通過梯度下降法求解二次方程來說明梯度下降法的效能。設求解的問題為 \({\bf{y}} = {\bf{Ax}}\)，那麼用梯度下降法求解 \({\bf{x}}^{*}\) 可以寫成優化問題： \({\bf{x}}^{*}=\mathop {\arg \min }\limits_{{\bf{x}} \in {\mathbb{R}^n}} \frac{1}{2}\left\| {{\bf{Ax}} - {\bf{y}}} \right\|_2^2\)。

4.1 Matlab code

• 梯度下降法函式

function [hat_x, error] = opt_gd(y,A,iter_max)
% gradient descent algorithm for solving y=Ax
% Louis Zhang 2020.12.10
% email:zhangyanfeng527@gmail.com
n = size(A,2) ;
x0 = zeros(n,1) ;
g0 = A'*(y-A*x0) ;
% 初始化梯度
% 初始化x_0
for k = 1:iter_max
    g1 = A'*(y-A*x0) ;
    % 計算梯度
    %  alpha = (g0'*g0)/(norm(A*g0)^2) ;
    alpha = 0.1 ;
    % 計算步長
    hat_x = x0 + alpha*g1 ;
    % 更新x_k+1
    error(k) = (norm(y-A*hat_x)/norm(y))^2 ;
    % error(k) = (norm(hat_x-x0)/norm(hat_x))^2 ;
    if error(k)<1e-8
        break;
    else
        x0 = hat_x ;
        g0 = g1 ;
    end
end
end

• 測試程式

clear
clc
N = 1000 ;
iter_max = 1000 ;
sig_max = 1 ;
% 最大奇異值
sig_min = 1 ;
% 最小奇異值
sig_num = linspace(sig_min,sig_max,N) ;
% 生成奇異值，最大奇異值為100，最小奇異值為1
V = diag(sig_num) ;
% 生成奇異值組成的對角矩陣
A_rand = randn(N,N) ;
% 生成隨機矩陣
A_orth = orth(A_rand) ;
% 隨機矩陣正交化
A = A_orth*V*A_orth^(-1) ;
% 得到設定條件數的矩陣A
cond(A)
% 試一下矩陣A的條件數是否正確
x = randn(N,1) ;
y = A*x ;

[hat_x,error] = opt_gd(y,A,iter_max) ;
% 使用梯度下降法求解y=Ax

figure(1)
plot(error,'m-o')
xlabel('迭代次數')
ylabel('NMSE')
title(['condition number=',num2str(sig_max/sig_min)])

4.2 模擬結果

圖3. 固定步長為0.5條件數為2時的結果 圖4. 最優步長，條件數為2時的結果 圖5. 最優步長，條件數為10時的結果

從結果可以看出來在相同條件數情況下，最優步長比固定步長收斂速度要快；在最優步長情況下，條件數越大收斂速度越慢。

5 討論

5.1 梯度下降法優點

• 梯度下降法的複雜度較低，比如在求解二次問題時，最小二乘的複雜度為 \(O\left( {{n^3}} \right)\)，而梯度下降法的複雜度為 \(O\left( {{n^2}} \right)\)。在求解大規模問題時優勢明顯。

5.2 梯度下降法缺點

• 梯度下降法的收斂速度較慢，原因是每次迭代時作為梯度方向作為尋優方向，梯度方向僅僅反映了點 \({\bf{x}}_k\) 的區域性性質，也就是說，對於區域性來說，梯度方向是最快的方向，而對於全域性來說，梯度方向不一定是最快的方向。
• 梯度下降法的收斂受初始點的影響較大，在非凸問題中，初始點如果選在最優點附近，則能以較快的速度收斂到全域性最優點，如果初始點與最優點較遠，則可能收斂到區域性最優點。
• 梯度下降法存在鋸齒收斂情況，尤其是在最優點附近時，由於梯度變化緩慢，收斂速度非常慢。

6 參考文獻

[1]]袁亞湘, 孫文瑜. 最優化理論與方法[M]. 科學出版社, 1997.
[2]https://www.cs.ubc.ca/~schmidtm/Courses/540-W18/L4.pdf
[3]https://www.cs.cmu.edu/~ggordon/10725-F12/slides/05-gd-revisited.pdf
[4]]http://users.ece.utexas.edu/~cmcaram/EE381V_2012F/Lecture_4_Scribe_Notes.final.pdf

更多精彩內容請關注訂閱號優化與演算法和加入QQ討論群1032493483獲取更多資料

往期精選：

• 經典正交多項式介紹及其應用

• SISO/MIMO通道估計（channel estimation）原理詳細圖解

• 毫米波大規模MIMO的AOA估計

• 線性降維：主成分分析PCA原理分析與模擬驗證

• 5G+AI：有哪些新的研究方向和新正規化？

• 簡述3D點雲配准演算法

• 5G為人工智慧與工業網際網路賦能|79頁高清PPT

• 智慧演算法|以動物命名的演算法

• 一份超全面的機器學習公共資料集

• 矩陣填充|奇異值閾值演算法

• 可重構/大規模智慧反射表面reconfigurable/large intelligent surface綜述

• 迭代硬閾值類演算法總結||IHT/NIHT/CGIHT/HTP

• 軟閾值迭代演算法（ISTA）和快速軟閾值迭代演算法（FISTA）

• 伍德伯裡矩陣恆等式（Woodbury matrix identity）

• 壓縮感知：一種新型亞取樣技術