最优化、凸优化、梯度下降和牛顿法

本文是对公众号SIGAI中凸优化，梯度下降，牛顿法等一系列文章的整理。

• SIGAI:理解凸优化
• SIGAI:理解梯度下降法
• SIGAI:理解牛顿法

最优化问题在机器学习中有非常重要的地位，很多机器学习算法最后都归结为求解最优化问题。凸优化（convex optimization）是最优化问题中非常重要的一类，对于机器学习来说，如果要优化的问题被证明是凸优化问题，则说明此问题可以被比较好的解决。

1. 最优化问题

最优化问题是求解函数极值的问题，包括函数的极大值和极小值（注意这里是极值，而不是最值，最值是全局的概念，如果一个函数有最值，那么只会有一个最值，而极值是从局部考虑的，一个还是的最值可能是极值也可能不是）。

比如求二次函数的最值，常使用配项法或者直接记住公式就可以了，而微积分为求解极值提供了统一的思路：函数的极值点处的导数为0。只要函数时可导的，就可以使用这个方法。

在机器学习中，我们一般将最优化问题统一表述为求解函数的极小值问题:$$mi{n}_{x}f(x)$$ 其中$x$为优化变量，$f$为目标函数。极大值问题可以转化为极小值问题求解
$$mi{n}_x-f(x)$$

1.1 可行域

通过为变量$x$添加约束条件，限定$x$的可行域，即满足约束条件的点构成的集合。

1.2 全局极小值（最小值）

一个优化问题的全局极小值是对可行域内所有的$x$，有：$$f(x^{\ast })\le f(x)$$

1.3 局部极小值

存在一个领域，对于在领域$||x-x^{\ast }||\le \delta$内的所有$x$，满足：$$f(x^{\ast })\le f(x)$$

求解一个一般性的最优化问题的全局极小值是非常困难的，至少要面临的问题是：函数可能有多个局部极值点，另外还有鞍点问题。对于第一个问题，我们找到了一个梯度为0的点，它是极值点，但不是全局极值，如果一个问题有多个局部极值，则我们要把所有局部极值找出来，然后比较，得到全局极值，这非常困难，而且计算成本相当高。第二个问题更严重，我们找到了梯度为0的点，但它连局部极值都不是，典型的是这个函数，在0点处，它的导数等于0，但这根本不是极值点：

2. 凸优化问题

我们常见的梯度下降法和牛顿法等基于导数的优化算法，转到的都是导数/梯度为0的点，而梯度等于0只是取得极值的必要条件而不是充分条件，如果要将这个必要条件变成充分条件，即：$$x点处的梯度等于0\Rightarrow x是问题的最优解$$ 问题将会得到简化。

如果对于问题加以限定：限定目标函数为凸函数；对于优化变量的可行域（注意，还要包括目标函数定义域的约束），我们限定它是凸集，这类问题有一个非常好性质，那就是局部最优解一定是全局最优解，我们称这类问题为凸优化问题。

2.1 凸集

2.1.1 凸集定义

对于$n$维空间中点的集合$C$，如果对集合中的任一两点$x$和$y$，以及实数$0\le \theta \le 1$，都有：$$\theta x+(1-\theta )y\in C$$则称该集合为凸集。相应的点$$\theta x+(1-\theta )y$$称为$x$和$y$的凸组合。下图中的左图为一个凸集，而右图不是。

2.1.1 凸集和约束条件的相关结论

1. 不带约束的优化问题，其优化变量的可行域是一个凸集

在$n$维向量空间$R^n$中，如果$x,y\in R^n$,显然有：$$\theta x+(1-\theta )y\in R^n$$

2. 一组线性等式约束条件，其确定的可行域是一个凸集

证明：给定$m$行$n$列的矩阵$A$和$m$维向量$b$, 仿射子空间定义为如下的向量空间集合 { x ∈ R n : A x = b } \{x \in R^n:Ax = b\} {x∈Rn:Ax=b} 假设$x_0,x_1\in R^n$为可行域中的两个解，则$$A{x}_0=b,A{x}_1=b$$ 对于任意的$0\le \theta \le 1$有：$$A(\theta x_0+(1-\theta )x_1)=\theta A{x}_0+(1-\theta )A{x}_1=\theta b+(1-\theta )b=b$$

3. 一组线性不等式约束条件，其确定的可行域是一个凸集

证明：给定$m$行$n$列的矩阵$A$和$m$维向量$b$, 仿射子空间定义为如下的向量空间集合 { x ∈ R n : A x ≤ b } \{x \in R^n:Ax \leq b\} {x∈Rn:Ax≤b} 假设$x_0,x_1\in R^n$为可行域中的两个解，则$$A{x}_0\le b,A{x}_1\le b$$ 对于任意的$0\le \theta \le 1$有：$$A(\theta x_0+(1-\theta )x_1)=\theta A{x}_0+(1-\theta )A{x}_1\le \theta b+(1-\theta )b=b$$

4. 多个凸集的交集还是凸集

2.2 凸函数

2.2.1 凸函数定义

在函数的定义域内，如果对于任意的$x$和$y$，以及实数$0\le \theta \le 1$，都满足如下条件：
$$f(\theta x+(1-\theta )y)\le \theta f(x)+(1-\theta )f(y)$$ 则函数为凸函数。这个不等式和凸集的定义类似。从图像上看，一个函数如果是凸函数，那么它是向下凸出去的。用直线连接函数上的任何两点A和B，线段AB上的点都在函数的上方，如下图所示：

如果上式中满足严格小于，则称函数为严格凸函数。

2.2.2 凸函数的判定

1. 一阶判定规则： $f(y)\ge f(x)+\nabla f(x{)}^T(y-x)$,为函数在任何点处的切线都位于函数的下方
2. 二阶判定规则：$f^{{}^{\prime \prime }}(x)\ge 0$,如果严格大于，则函数时严格凸的。
   对于多元函数，如果它是凸函数， 则其Hessian矩阵为半正定矩阵。如果Hessian矩阵是正定的，则函数是严格凸函数。

根据多元函数极值判别法，假设多元函数在点M的梯度为0，即M是函数的驻点，则有：

• 如果Hessian矩阵正定，函数在该点有极小值
• 如果Hessian矩阵负定，函数在该点有极大值
• 如果Hessian矩阵不定，还需要看更高阶的导数

从单变量的角度来说，二阶导数一直大于0表示一阶导数的值一直在增加，又有一阶导数有0值，即一阶导数必然是由负值增加到0，在由0增加到正值，那么原函数必然是先降低在上升，函数必有极小值。

对于n阶矩阵$A$，对于任意的非0的n维向量$x$, 都有$x^{T}Ax>0$，则称矩阵A（A可以看做二次型的系数）为正定矩阵。判定矩阵正定的常用方法有以下几种：

• 矩阵的特征值全大于0
• 矩阵的所有顺序主子式都大于0
• 矩阵合同于单位阵$I$

$x^{T}Ax<0$则称矩阵为负定矩阵，$x^{T}Ax\ge 0$则矩阵为半正定。

凸函数的非负线性组合时凸函数，假设$f_i$是凸函数，并且$w_i\ge 0$，则：$f(x)={\sum }_{i=1}^k{w}_i{f}_i(x)$是凸函数
注意：这边是函数的非负性组合，不是特征的非负性组合

2.3 凸优化

如果一个最优化问题的可行域是凸集，并且目标函数是凸函数，则该问题为凸优化问题。凸优化问题可以形式化的写成：$$minf(x)$$$$x\in C$$ 其中$x$为优化变量，$f$为凸目标函数，$C$为优化变量的可行域，为一个凸集。

将可行域具体化，凸优化问题可以表述为：$$minf(x)$$$$g_i(x)\le 0,i=1,...,m$$$$h_i(x)=0,i=1,...,p$$
凸优化问题有一个重要的特性：所有局部最优解都是全局最优解。这个特性可以保证我们在求解时不会陷入局部最优解，即如果找到了问题的一个局部最优解，则它一定也是全局最优解，这极大的简化了问题的求解。

3. 最优化问题求解

从上文的介绍中，我们知道函数的极值点处，一阶导数为0，根据二阶导数的Hessian矩阵，可以判定该值是极小值，极大值或者是鞍点，那么对于最优化问题，直接求导，令梯度为0，然后解方程不就可以了吗？事实上没有那么简单，看下面这个复杂函数：
$$f(x,y)=x^3-2X^2+e^{xy}-y^3+10y^2+100sin(xy)$$ 分别对x,y求导，并令其为0时，得到：$$\{\begin{array}{ll}3x^2-4x+y{e}^{xy}+100ycos(xy)=0& \\ x{y}^{xy}-3y^2+20y+xcos(xy)=0& \end{array}$$ 这个方程非常难以求解，对于有指数函数，对数函数，三角函数的方程，我们称为超越方程，求解的难度并不比求极值本身小。

3.1 迭代法求解思想

精确的求解不太可能，因此只能求近似解，这称为数值计算。工程上实现时通常采用的是迭代法，它从一个初始点开始，反复使用某种规则从移动到下一个点，构造这样一个数列，直到收敛到梯度为0的点处。即有下面的极限成立$$li{m}_{k\to +\infty }\nabla f(x_k)=0$$其中k为迭代次数。
这些规则一般会利用一阶导数信息即梯度；或者二阶导数信息即Hessian矩阵。这样迭代法的核心是得到这样的由上一个点确定下一个点的迭代公式：$$x_{k+1}=h(x_k)$$

3.2 梯度下降法

梯度下降法只利用了函数的一阶导数信息。

3.2.1 梯度下降法原理

如果一元函数$f(x)$n阶可导，它的泰勒展开公式为：$$f(x+\Delta x)=f(x)+f^{\prime }(x)\Delta x+\frac{1}{2}{f}^{\prime \prime }(x)(\Delta x{)}^2+...+\frac{1}{n!}{f}^{(n)}(x)(\Delta x{)}^n$$ 多元函数的泰勒展开为：$$f(X+\Delta X)=f(X)+(\nabla f(X){)}^T\Delta x+o(\Delta X)$$ 此处忽略了二次及更高项。$$f(X+\Delta X)-f(X)=(\nabla f(X){)}^T\Delta x+o(\Delta X)$$ 如果$\Delta X$足够小，$o(\Delta X)$可以忽略：$$f(X+\Delta X)-f(X)\approx (\nabla f(X){)}^T\Delta X$$
因为$X$为多元变量，不像一元变量$x$只有左右两个方向，这是$X$可以走的方向有很多，该怎么走呢，因为我们想要$f(X)$朝极小值走，即$f(x+\Delta X)<f(X)$，所以只要保证$(\nabla f(X){)}^T\Delta X<0$即可。因为$$(\Delta f(X){)}^T\Delta X=\Vert \nabla f(X)\Vert \Vert \Delta X\Vert cos\theta$$ 其中$\Vert \cdot \Vert$表示向量的模，$\theta$为向量之间的夹角。因为向量的模一定大于等于0，所以只需要$cos\theta \le 0$就能得到我们想要的结果。

因为$cos\theta =[-1,1]$，当$\theta =\pi$时有极小值-1，此时$f(X+\Delta X)$沿着$f(X)$的梯度方向，达到了最大的下降值，该下降的值为：$$\Vert \nabla f(X)\Vert \Vert \Delta X\Vert$$
设$\Delta X=-\alpha \nabla f(X)$，$\alpha$是一个接近于0的整数，称为步长，由人工设定，用于保证$x+$在$x$的邻域，此时有：$$(\nabla f(X){)}^T\Delta X=-\alpha (\nabla f(X){)}^T(\nabla f(X))\le 0$$ 而变量$x$的迭代公式为$$X_{k+1}=X_k-\alpha \nabla f(X_k)$$

只要没有到达梯度为0的点，则函数值会沿着负梯度递减，最终会收敛到梯度为0的点，这就是梯度下降法。迭代终止的条件是函数的梯度值为0（实际实现时是接近于0），此时认为已经达到极值点。注意我们找到的是梯度为0的点，这不一定就是极值点，后面会说明。梯度下降法只需要计算函数在某些点处的梯度，实现简单，计算量小。

3.2.2 梯度下降法实现细节

1. 初始值的设定

一般的，对于不带约束条件的优化问题，我们可以将初始值设置为0，或者设置为随机数，对于神经网络的训练，一般设置为随机数，这对算法的收敛至关重要

2. 学习率的设定

学习率设置为多少，也是实现时需要考虑的问题。最简单的，我们可以将学习率设置为一个很小的正数，如0.001。另外，可以采用更复杂的策略，在迭代的过程中动态的调整学习率的值。比如前1万次迭代为0.001，接下来1万次迭代时设置为0.0001

3. 梯度下降的问题

• 局部极小值
• 鞍点

3.2.3 梯度下降法的变种

梯度下降法有大量的变种，它们都只利用之前迭代时的梯度信息来构造每次的更新值

随机梯度下降

每次使用一个样本对参数进行更新，单个样本的梯度方向不一定是最优的，所以和梯度下降相比，随机梯度下降需要更多的迭代次数达到收敛效果。
但是随机梯度下降因为下降方向不定，所以可能会跳出鞍点和当前邻域内局部最优值。

批量梯度下降

批量梯度下降结合了梯度下降和随机梯度下降，每次使用一个小批量的训练数据对参数进行迭代，特别是对不能完全load到内存的训练数据，建议使用批量梯度下降。

动量

动量累计了之前的权重更新值，假设$V_{t+1}$为$t+1$次迭代的动量，优化变量$x$的迭代方程为：$$X_{t+1}=X_t+V_{t+1}$$ $V_{t+1}$取代了之前的梯度项。动量项的计算公式为：$$V_{t+1}=-\alpha \nabla f(x_t)+\mu V_t$$动量项累积了之前的梯度信息，类似于保持行走时的惯性，以避免来回震荡，加快收敛速度。

AdaGrad

AdaGrad为自适应梯度，即adaptive gradient算法，AdaGrad根据前几轮迭代时的历史梯度值来调整学习率，参数更新公式为：$$(x_{t+1}{)}_i=(x_t{)}_i-\frac{\alpha }{\sqrt{{\sum }_{i=1}^t(f_j^{\prime }(x_i){)}^2+ϵ}}{f}_t^{\prime }(x_i)$$ $\alpha$为全局学习步长，$f_t(x_i)$是第$i$个优化变量的第$t$次迭代梯度，$ϵ$是一个很小的整数，避免出现分母为0的情况。

每个变量都随着学习的进行，根据历史学习率累积总量来决定当前学习率减小的程度。在机器学习的应用中，AdaGrad非常适合样本稀疏的问题，因为稀疏的样本下，每次梯度下降的方向，以及涉及的变量都可能有很大的差异。AdaGrad的缺点是虽然不同变量有了各自的学习率，但是初始的全局学习率还是需要手工指定。如果全局学习率过大，优化同样不稳定；而如果全局学习率过小，因为AdaGrad的特性，随着优化的进行，学习率会越来越小，很可能还没有到极值就停滞不前了。

Adam

Adam算法全称为adaptive moment estimation。它用梯度项构造了两个向量$m$和$v$，这两个向量的更新公式为$$(m_t{)}_i={\beta }_1(m_{t-1}{)}_i+(1-{\beta }_1)(f_t^{\prime }{)}_i$$$$(v_t{)}_i={\beta }_2(v_{t-1}{)}_i+(1-{\beta }_2)(f_t^{\prime }{)}_i^2$$ 其中${\beta }_1,{\beta }_2$由人工指定，$i$为优化向量的分量下标，参数的更新公式：
$$(x_{t+1}{)}_i=(x_t{)}_i-\alpha \frac{\sqrt{1-({\beta }_2{)}_i^2}}{1-({\beta }_1{)}_i^t}\frac{(m_t{)}_i}{\sqrt{(v_t{)}_i}+ϵ}$$
这里，$m$代替梯度，$v$构造学习率。

除了以上介绍的随机梯度改进方法，还有如AdaDelta、 NAG等其他的一些方法。

3.3 牛顿法

第三部分使用了一阶泰勒展开式来解释梯度下降的原理，并使用了一阶导数迭代求参。牛顿法则使用到了二阶导数近似$x$邻域内的值。

3.3.1 牛顿法推导

多元函数在$X_0$的二阶泰勒展开为：
$$f(X)=f(X_0)+(\nabla f(X_0){)}^T(X-X_0)+\frac{1}{2}(X-X_0{)}^T({\nabla }^{2}f(X_0))(X-X_0)+o((X-X_0{)}^2)$$
忽略二次及以上项。并对上式两边同时求梯度：
$$\nabla (X)=\nabla f(X_0)+{\nabla }^{2}f(X_0)(X-X_0)$$
其中${\nabla }^{2}f(X_0)$为Hessian矩阵，记为$H$。 令函数的梯度为0，则有：
$$X=X_0-({\nabla }^{2}f(X_0){)}^{-1}\nabla f(X_0)$$
将梯度向量简写为$g$, 上面的公式简写为：
$$X=X_0-H^{-1}g$$
从初始点$X_0$处开始，反复计算函数的Hessian矩阵和梯度向量，然后用下式进行迭代：
$$X_{k+1}=X_k-H_k^{-1}{g}_k$$
其中$-H^{-1}g$称为牛顿方向，迭代终止的条件是梯度的模接近于0，或者函数值下降小于制定阈值。

3.3.2 实现细节

其中$\gamma$是一个超参，接近于0，该参数保证了$X_{k+1}$在$X_k$的邻域内，从而可以忽略泰勒展开的高次项。如果目标函数是二次函数，Hessian矩阵是一个常数矩阵，对于任意给定的初始点，牛顿法只需要一步迭代就可以收敛到极值点。

以$x^2+y^2$为例：

如果我们把步长设置的足够大（在这里为1），则算法一步就收敛了。在这里，初始迭代位置为(0,4)，最优解为(0,0)。

对于不带约束条件的问题，我们可以将x的初始值设定为任意值，最简单的，可以设置为全0的向量。迭代终止的判定规则和梯度下降法相同，是检查梯度是否接近于0。

和梯度下降法相比牛顿法有更快的收敛速度，但每一步迭代的成本也更高。在每次迭代中，除了要计算梯度向量还要计算Hessian矩阵，并求解Hessian矩阵的逆矩阵。实际实现时一般不直接求Hessian矩阵的逆矩阵，而是求解如下方程组：
$H_{k}d=-g_k$
求解这个线性方程组一般使用迭代法，如共轭梯度法等。

3.3.3 存在的问题

和梯度下降法一样，牛顿法寻找的也是导数为0的点，这不一定是极值点，因此会面临局部极小值和鞍点问题，这在之前的文章中已经介绍过了，这里不再重复。

牛顿法面临的另外一个问题是Hessian矩阵可能不可逆，从而导致这种方法失效。此外，求解Hessian矩阵的逆矩阵或者求解线性方程组计算量大，需要耗费大量的时间。为此，提出了拟牛顿法这种改进方案，在后面会介绍。

除此之外，牛顿法在每次迭代时序列xi可能不会收敛到一个最优解，它甚至不能保证函数值会按照这个序列递减。解决第一个问题可以通过调整牛顿方向的步长来实现，目前常用的方法有两种：直线搜索和可信区域法。可信域牛顿法在后面也会介绍。

直线搜索
牛顿法并不能保证每一步迭代时函数值下降，也不保证一定收敛。为此，提出了一些补救措施，其中的一种是直线搜索（line search）技术，即搜索最优步长。具体做法是让取一些典型的离散值，如0.0001,0.001,0.01等，比较取哪个值时函数值下降最快，作为最优步长。

可信域牛顿法(不是很理解)
可信域牛顿法（Trust Region Newton Methods）可以求解带界限约束的最优化问题，是对牛顿法的改进。在可信域牛顿法的每一步迭代中，有一个迭代序列$x_k$，一个可信域的大小${\nabla }_k$，以及一个二次目标函数：
$q_k(s)=(\nabla f(X_k){)}^{T}s+\frac{1}{2}{s}^T{\nabla }^{2}f(x_k)s$
这个式子可以通过泰勒展开式得到，忽略二次以上的项，这是对函数下降值：$f(x_k+s)-f(x_k)$的近似。算法寻找一个$s_k$，在满足约束条件$||s||\le {\Delta }_k$下近似最小化$q_k(s)$。接下来根据以下式子来更新$x_k$和${\Delta }_k$:
$$p_k=\frac{f(x_k+s_k)-f(x_k)}{q_k(s_k)}$$
这是函数值的实际减少量和二次近似模型预测方向导致的函数减少量的比值。迭代方向可以接受的条件是$p_k$足够大，由此可以得到参数的更新规则为：
$$x_{k+1}=\{\begin{array}{llcc}{x}_k+s_k,& & p_k>{\eta }_0& \\ x_k,& & p_k\le {\eta }_0& \end{array}$$
其中${\eta }_0$是一个人工设定的值。${\Delta }_k$的更新规则取决于人工设定的正常数${\eta }_1$和${\eta }_2$,其中：${\eta }_1<{\eta }_2<1$.

拟牛顿法
牛顿法在每次迭代时需要计算出Hessian矩阵，然后求解一个以该矩阵为系数矩阵的线性方程组，这非常耗时，另外Hessian矩阵可能不可逆。为此提出了一些改进的方法，典型的代表是拟牛顿法（Quasi-Newton）。

拟牛顿法的思想是不计算目标函数的Hessian矩阵然后求逆矩阵，而是通过其他手段得到Hessian矩阵或其逆矩阵的近似矩阵。具体做法是构造一个近似Hessian矩阵或其逆矩阵的正定对称矩阵，用该矩阵进行牛顿法的迭代。将函数在点处进行泰勒展开，忽略二次以上的项，有：

这个条件称为拟牛顿条件，用来近似代替Hessian矩阵的矩阵需要满足此条件。根据此条件，构造出了多种拟牛顿法，典型的有DFP算法、BFGS算法、L-BFGS算法等，在这里我们重点介绍BFGS算法。下图列出了常用的拟牛顿法的迭代公式（图片来自维基百科）：

BFGS算法是它的四个发明人Broyden，Fletcher，Goldfarb和Shanno名字首字母的简写。算法的思想是构造Hessian矩阵的近似矩阵：

总结：
梯度下降，牛顿法等可以解所有最优化问题，但是凸优化问题得到的解是全局最优的，而其他的最优化问题得到的解则不一定是全局最优。