梯度下降法原理与仿真分析||系列（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获取更多资料

往期精选：