引言

本节将介绍无约束优化问题的常用求解方法，包括坐标轴交替下降法、最速下降法。

本节是对优化算法(十~十七)最速下降法(梯度下降法)的理论补充，其中可能出现一些定理的证明过程这里不再赘述，并在相应位置附加链接。

总体介绍

从本节开始，将介绍四大类无约束优化问题的常用求解方法：

坐标轴交替下降法；
最速下降法；
牛顿法；
拟牛顿法。

这些方法的核心区别在于：下降方向选择策略的差异性。通过介绍各算法选择下降方向的方式，并延伸至该算法的特点。

回顾：线搜索下降算法

关于最小化目标函数 $\min f(x)$ 的无约束优化问题，线搜索下降算法的迭代步骤表示如下：

给定数值解序列 $\{x_k\}_{k=0}^{\infty}$ 的迭代初始点 $x_0(k=0)$ ；
这仅是从数学角度对数值解序列进行描述。如果从算法角度，它不可能是一个长度为无穷大的序列。可以通过终止条件使迭代算法停止。
判断点 $x_k$ 是否满足终止条件：是，则终止；
寻找 $x_k$ 位置的下降方向 $\mathcal D_k$ ；
选择合适的步长 $\alpha_k \geq 0$ ，使得：
$f(x_k + \alpha_k \cdot \mathcal D_k) < f(x_k)$
令： $x_{k+1} = x_k + \alpha_k \cdot \mathcal D_k$ ；并令 $k = k + 1$ ，转步骤 $2$ 。

其中：

常用终止条件： $\|\nabla f(x_k)\| \leq \epsilon$
- 其中 $\epsilon$ 是一个较小的正值。例如 $10^{-6}$ 。
- 如果满足该条件，意味着： $x_k$ 点处的梯度 $\nabla f(x_k)$ 已经充分接近于 $0$ 。
步长选择方式：基于区间的直接搜索法；非精确搜索准则(五~七)；
包括 $\text{Armijo,Glodstein,Wolfe}$ 准则。因为仅仅让 $\{f(x_k)\}_{k=0}^{\infty}$ 收敛并不是其达到最优解的充要条件。详见线搜索方法(步长角度；非精确搜索)
下降方向；

针对不同的下降方向选择方式，产生不同种类的算法。而我们更关心的是对应算法产生的数值解序列 $\{x_k\}_{k=0}^{\infty}$ 是否能够收敛至最优解 $x^*$ ，如果能够收敛至最优解 $x^*$ ，需要关心它的收敛速度情况。

收敛速度的衡量方式

对应文章详见：优化算法(九)收敛速度的简单认识

线性收敛范围

假设数值解序列 $\{x_k\}_{k=0}^{\infty} \Rightarrow x^*$ ，如果存在极限：
很明显，关于 $\beta$ 的取值范围： $\beta \in [0,1]$ 。

其中当 $\beta=1$ 时，分母与分子之间的差异性可视作完全相同;换句话说，当 $k$ 充分大时,两者之间的差距确实存在，但小到可以忽略不计。称这种收敛方式为次线性收敛。
当 $0<\beta<1$ 时，可以明显观察到分母与分子之间存在比值的大小关系；通过该比值 $\beta$ 可以明显观察到迭代过程中呈线性的收敛效果。
当 $\beta = 0$ 时，和 $\beta = 1$ 相反，当 $k$ 充分大时，分母与分子之间的差距足够大,甚至分子与分母相比，小到可以忽略不计。
$\mathop{\lim}\limits_{k \Rightarrow \infty} \frac{\|x_{k+1} - x^*\|}{\|x_k - x^*\|} = \beta$

根据 $\beta$ 的不同取值，有：

当 $\beta < 1$ 时，称数值解序列 ${x_k\}$ 为线性收敛；
当 $\beta = 0$ 时，则称数值解序列 ${x_k\}$ 为超线性收敛。

示例：假设 $\begin{aligned}\beta = \frac{1}{2}\end{aligned}$ ，那么：
$\begin{cases} \begin{aligned} \|x_{k+1} -x^*\| & = \frac{1}{2} \|x_k - x^*\| \\ \|x_{k+2} - x^*\| & = \frac{1}{2} \|x_{k+1} - x^*\| = \frac{1}{4}\|x_k - x^*\| \\ \vdots \\ \end{aligned} \end{cases}$
可以明显观察到其呈线性的收敛效果。

高阶收敛范围

如果存在 $\geq 1$ ，有：
$\mathop{\lim}\limits_{k \Rightarrow \infty} \frac{\|x_{k+1} - x^*\|}{\|x_k - x^*\|^p} = \beta < +\infty$
则称 ${x_k\}$ 为 $p$ 阶收敛。

牛顿法在适当条件下被证明是二阶收敛。
可以想象，当 $p > 1$ 时，相比于线性收敛范围,高阶收敛必然是更高级别的收敛速度。从而有如下表达：
当 $p > 1$ 时， $p$ 阶收敛必然为超线性收敛，但反之不一定成立。

验证：当 $p > 1$ 时，可以将上式拆解为如下形式：
$\mathop{\lim}\limits_{k \Rightarrow \infty} \frac{\|x_{k+1} - x^*\|}{\|x_k - x^*\|^p} = \mathop{\lim}\limits_{k \Rightarrow \infty} \left(\frac{\|x_{k+1} - x^*\|}{\|x_k - x^*\|} \cdot \frac{1}{\|x_k - x^*\|^{p-1}}\right)$

其中第一项描述的是线性收敛范围；观察第二项： $\begin{aligned}\lim_{k \Rightarrow \infty} \frac{1}{\|x_k - x^*\|^{p-1}}\end{aligned}$ 在 $p > 1$ 条件下，其结果是 $+\infty$ 。
如果需要 $\begin{aligned}\mathop{\lim}\limits_{k \Rightarrow \infty} \frac{\|x_{k+1} - x^*\|}{\|x_k - x^*\|} \cdot \infty = \beta < \infty\end{aligned}$ ，必然需要 $\begin{aligned}\mathop{\lim}\limits_{k \Rightarrow \infty} \frac{\|x_{k+1} - x^*\|}{\|x_k - x^*\|} = 0\end{aligned}$ ，即超线性收敛。

二次终止性

关于判断一个算法的优劣性，除去收敛速度这个评价标准外，优化问题本身也可以作为算法优劣性的评价标准。算法针对某类简单问题：

可能无法在有限迭代步骤内实现收敛；
可能会在有限迭代步骤内实现收敛，但计算代价过大；

这样的算法本身存在问题。相反，如何衡量简单问题的基准 $?$ 通常将目标函数为凸二次函数作为基准：
矩阵 $\mathcal Q$ 至少是半正定矩阵。
$\frac{1}{2}x^T \mathcal Qx + \mathcal C^T x \quad \mathcal Q \succcurlyeq 0$
如果针对上述问题在有限迭代步骤内接近最优解，我们称该算法具有二次终止性。

朴素算法：坐标轴交替下降法

其基本思想表示为：给定初始点 $x_0 \in \mathbb R^n$ ，依次沿坐标轴 $e_1,e_2,\cdots,e_n$ 进行搜素。

关于坐标轴交替下降法,它并不想在迭代步骤中花费代价计算下降方向,而是直接选择坐标轴方向作为下降方向。
这与吉布斯采样方法的思想——坐标上升法如出一辙。

对应算法框架表示如下：

给定初始点 $x_0;k=0;$
依然判断 $\|\nabla f(x_k)\| \leq \epsilon$ ：如果满足，终止；
记 $y_0 = x_k$ ，令：
$\begin{cases} y_i = y_{i-1} + \alpha_i \cdot e_i \\ \alpha_i = \mathop{\arg\min} f(y_{i-1} + \alpha \cdot e_i) \quad i=1,2,\cdots,n \end{cases}$
解释：实际上该步骤是一个 $n$ 次循环。这里的 $y_i(i=1,2,\cdots,n)$ 分别表示特征空间中的具体点。这里以二维特征 $x_k \in \mathbb R^2 \Rightarrow (e_1,e_2)$ 为例，使用图像描述该过程：
- 初始状态下， $y_0 = x_k:(x_1^{(k)},x_2^{(k)})$ ；
- 在除去 $e_1$ 外，其他维度固定的条件下，此时固定优化方向为 $e_1$ ，在该方向上的最优步长 $\alpha_1$ 可表示为关于步长变量 $\alpha$ 函数 $\phi(\alpha)$ 的最优解：
  $\alpha_1 = \mathop{\arg\min}\limits_{\alpha} \phi(\alpha) = \mathop{\arg\min}\limits_{\alpha} f(y_0 + \alpha \cdot e_1)$
- 找到 $\alpha_1$ 后，通过 $y_1 = y_0 + \alpha_1 \cdot e_1$ 可以得到第一次循环结束后更新的位置；
- 同上，继续循环，寻找除去 $e_2$ 外，其他维度固定的条件下，求出 $e_2$ 方向上的最优步长 $\alpha_2$ ，以此类推。直到 $n$ 个维度全部被遍历一次为止，得到 $y_n= x_{k+1}$ 。对应图像表示如下：
  当然这里 $n = 2$ 。
在得到 $x_{k+1} = y_n$ 后， $k = k + 1$ ，并步骤 $2$ ，直到满足条件为止。

该算法的优势在于：

不需要花费额外代价计算下降方向；
步骤 $3$ 的循环中， $e_i \in \mathbb R(i=1,2,\cdots,n)$ ，因而计算上相对简单。
当目标函数 $f (x)$ 中的决策变量 $\in \mathbb R^n$ ，其各分量 $x_i(i=1,2,\cdots,n)$ 之间的交叉程度很小时，该算法框架会非常有效。
什么是交叉程度很小——可理解为各分量之间的关联关系较小，甚至是线性无关。例如各分量满足可分离函数：各分量各算各的~
$\begin{aligned} \min f(x) & = \min [f_1(x_1) + f_2(x_2)+\cdots + f_n(x_n)] \\ & = \sum_{i=1}^n \min f_1(x_1) \end{aligned}$

相反，该算法的劣势在于：对于一般问题，该算法得到的数值解序列 $\{x_k\}_{k=0}^{\infty}$ 不一定收敛。
如果决策变量内各分量之间的关联性程度较高，其产生的结果并不容易收敛，吉布斯采样同样存在这种缺陷。

一种改进方法描述：将线搜索方法与坐标轴交替下降法交替使用从而使数值解序列收敛。具体改进步骤如下：
前面步骤并没有发生变化，在通过坐标轴交替下降法找到 $\bar{x}_k$ 后，能够确定： $f(\bar{x}_k) \leq f(x_k)$ ,也就是说： $x_k \Rightarrow \bar{x}_k$ 的方向 $\mathcal D_k$ 一定是下降方向。

给定初始点 $x_0;k=0;$
依然判断 $\|\nabla f(x_k)\| \leq \epsilon$ ：如果满足，终止；
记 $y_0 = x_k$ ，令：
$\begin{cases} y_i = y_{i-1} + \alpha_i \cdot e_i \\ \alpha_i = \mathop{\arg\min} f(y_{i-1} + \alpha \cdot e_i) \quad i=1,2,\cdots,n \end{cases}$
从而得到 $\bar{x}_{k}$ 。
以 $\bar{x}_k$ 为起始点， $\mathcal D_k:x_k \Rightarrow \bar{x}_k$ 为下降方向使用线搜索方法选择合适步长，从而得到新的更新结果 $x_{k+1}$ ；
依然是基于 $2$ 维特征，对应示例图像表示如下。
得到 $x_{k+1}$ 后， $k = k + 1$ ，并返回步骤 $2$ 。

最速下降法(梯度下降法)

其基本思想表示为：在迭代过程中，选择 $x_k$ 处的负梯度方向作为搜索方向。即： $\mathcal D_k = - \nabla f(x_k)$ 。
而负梯度方向也被称作最速下降方向。

从泰勒展开式的角度观察，根据线搜索方法(方向角度)的下降方向的推导过程可知：若判断 $x_k$ 处的某方向 $\mathcal D$ 是否为下降方向，只需判断：
$[\nabla f(x_k)]^T \mathcal D < 0$
那么方向 $\mathcal D$ 就是 $x_k$ 位置的下降方向。当 $\mathcal D = -\nabla f(x_k)$ 时，能够使 $[\nabla f(x_k)]^T \mathcal D$ 达到最小值：
这里仅关注向量 $\nabla f(x_k),\mathcal D$ 的方向信息，因而设 $\|\nabla f(x_k)\| = \|\mathcal D\| = 1$ 。
$[\nabla f(x_k)]^T \mathcal D = \|\nabla f(x_k)\| \cdot \|\mathcal D\| \cos \theta$
其中 $\theta$ 表示向量 $\nabla f(x_k),\mathcal D$ (不分先后)之间的夹角。当 $\mathcal D,\nabla f(x_k)$ 之间夹角为 $\begin{aligned}\frac{\pi}{2}\end{aligned}$ 时，能够取到 $\cos \theta$ 的最小值 $- 1$ 。
如果从方向导数的角度观察： $[\nabla f(x_k)]^T \mathcal D$ ，它可以看作： $x_k$ 所在位置处关于 $\mathcal D$ 的方向导数。在凸函数铺垫：梯度与方向导数中介绍过，对应方向导数可表示为：
这里示例 $x_k$ 是二维特征，坐标为 $(x, y)$ 。
$\begin{aligned} \frac{\partial \mathcal Z}{\partial \mathcal D}\mid_{(x,y)} & = f_x(x_k) \cdot \cos \alpha + f_y(x_k) \cdot \cos \beta \\ & = \underbrace{[f_x(x_k),f_y(x_k)]}_{[\nabla f(x_k)]^T} \begin{pmatrix} \cos \alpha \\ \cos \beta \end{pmatrix} \\ & = [\nabla f(x_k)]^T \mathcal D \end{aligned}$
关于方向导数的性质：
这意味着： $[\nabla f(x_k)]^T \mathcal D$ 达到最小值，意味着函数值下降的越剧烈。
- 当 $[\nabla f(x_k)]^T \mathcal D > 0 \Rightarrow$ 在 $x_k$ 位置沿着 $\mathcal D$ 方向的函数值上升；反之， $[\nabla f(x_k)]^T \mathcal D < 0 \Rightarrow$ 在 $x_k$ 位置沿着 $\mathcal D$ 方向的函数值下降。
- $|\nabla f(x_k)^T \mathcal D|$ 越大 $\Rightarrow$ 上升/下降的越猛烈；反之， $|\nabla f(x_k)^T \mathcal D|$ 越小 $\Rightarrow$ 上升/下降的越平缓。

梯度下降法的特点

优点：
梯度下降法能够收敛，并且其下降方向被指定为负梯度方向 $-\nabla f(x_k)$ 。

缺陷：

收敛速度慢，即便是在凸函数甚至是强凸函数最快也只能达到线性收敛；
相关证明见：梯度下降法在强凸函数上的收敛性证明以及梯度下降法在凸函数上的收敛性。

归纳：
- 梯度下降法仅使用负梯度方向作为搜索方向，换句话说：在考虑搜索方向的过程中，仅考虑了一阶梯度 $\nabla f(\cdot)$ 信息；实际上，二阶梯度信息 $(\text{Hessian Matrix})$ 也可以用来判断搜索方向。
- 其次，假设在最速下降法的过程中，由于方向 $\mathcal D_k$ 已被确定，那么最优步长 $\alpha_k$ 是关于 $\phi(\alpha) = f(x_k + \alpha \cdot \mathcal D_k)$ 的精确最小点：
  $\alpha_k = \mathop{\arg\min}\limits_{\alpha} \phi(\alpha) =\mathop{\arg\min}\limits_{\alpha} f(x_k + \alpha \cdot \mathcal D_k)$
  令 $\phi'(\alpha) \triangleq 0$ ，必然有：
  $\phi'(\alpha_k) = [\nabla f(x_k + \alpha_k \cdot \mathcal D_k)]^T \mathcal D_k = [\nabla f(x_{k+1})]^T[-\nabla f(x_k)] = 0$
  这意味着：梯度向量 $\nabla f(x_{k+1})$ 与梯度向量 $\nabla f(x_k)$ 垂直。
  而这个垂直于 $\mathcal Z$ 字形的缺陷是同一个缺陷：它仅能在迭代步骤中找到局部最优方向,而不是全局最优方向。也就是说：梯度下降法是一个贪心算法。
$\text{ZigZag}$ 现象：在迭代过程中，其收敛路径呈 $\mathcal Z$ 字形；
见下方代码示例与图像。可以看出：其搜索路径呈线 $\mathcal Z$ 字形，并且每一次迭代的方向均不是全局最优。
不具备二次终止性，也就是说：关于凸二次函数的最优化问题，仅仅通过有限次迭代步骤，无法收敛至最优解。

针对最速下降法缺陷代码示例

针对梯度下降法上述缺陷问题，以凸二次函数的最优化问题： $\begin{aligned}\min f(x,y) = \frac{1}{2} x^2 + 2 y^2\end{aligned}$ 为例，使用最速下降法近似求解最优解。对应代码表示如下：

import numpy as np
import math
import matplotlib.pyplot as plt


def f(x,y):
    return 0.5 * (x ** 2) + 2 * (y ** 2)

def ConTourFunction(x,Contour):
    return math.sqrt(0.5 * (Contour - (0.5 * (x ** 2))))

def Derfx(x):
    return x

def Derfy(y):
    return 4 * y

def GradientDescent(stepTime=10,epsilon=0.1):
    Start = (2.0,1.0)
    LocList = list()
    LocList.append(Start)

    for _ in range(stepTime):

        DerStart = (Derfx(Start[0]),Derfy(Start[1]))
        for step in list(np.linspace(0.0,1.0,1000)):
            Next = (Start[0] - (DerStart[0] * step),Start[1] - (DerStart[1] * step))
            DerfNext = Derfx(Next[0]) * (-1 * DerStart[0]) + Derfy(Next[1]) * (-1 * DerStart[1])

            if abs(DerfNext) <= epsilon:
                LocList.append(Next)
                Start = Next
                epsilon /= 5.0
                break

    ContourList = [0.1,0.2,0.5,1.0]
    LimitParameter = 0.0001
    plt.figure(figsize=(10,5))
    for Contour in ContourList:
    	# 设置范围时，需要满足x的定义域描述。
        x = np.linspace(-1 * math.sqrt(2 * Contour) + LimitParameter,math.sqrt(2 * Contour) - LimitParameter,200)
        y1 = [ConTourFunction(i,Contour) for i in x]
        y2 = [-1 * j for j in y1]
        plt.plot(x,y1,'--',c="tab:blue")
        plt.plot(x,y2,'--',c="tab:blue")

    plotList = list()
    for (x,y) in LocList:
        plotList.append((x,y))
        plt.scatter(x,y,s=50,facecolor="none",edgecolors="tab:red",marker='o')
        if len(plotList) < 2:
            continue
        else:
            plt.plot([plotList[0][0],plotList[1][0]],[plotList[0][1],plotList[1][1]],c="tab:red")
            plotList.pop(0)
            
    plt.plot([0,2],[0,1],'--',c="tab:green")
    plt.show()

if __name__ == '__main__':
    GradientDescent()