梯度下降法在 OI 中的应用

· · 算法·理论

梯度下降法（Gradient descent）是一种用于求解可微函数局部极小值的迭代算法，常用于机器学习中 [^1]。

[^1]: Gradient descent - Wikipedia

下面，我们将简要介绍梯度下降的过程，并梳理其在 OI 中的应用。若有不当之处，恳请指出，不胜感激！

若神犇发现了更好的应用，请告知我，我随时更新本文！！

前置知识

相信各位已经对函数和向量的相关概念有基本的了解。不过，要弄清楚梯度下降法是怎么一回事，还有几个核心概念需要了解。在这里，我尽量用最简单易懂的语言介绍它们。若您已经熟悉这些概念，可跳过此部分。

最优化

最优化 (Mathematical optimization) 指从一组可选择的方案中，根据一定标准选择最佳方案的过程。一般分为离散优化、连续优化两个子领域 [^2]。

[^2]: Mathematical optimization - Wikipedia

在梯度下降中，我们主要研究连续优化 (Continuous Optimization) 问题：求解连续函数的极小值。这个函数称为目标函数 (Objective function) 。

在 OI 中有很多离散优化问题，一般是从有限的自变量取值中找到某个函数的极小/极大值。梯度下降法不负责解决此类问题，建议考虑模拟退火。

无约束优化

我们遇到的连续优化问题都有必然的约束，即自变量必须在函数定义域内（大多数情况下都是 \boldsymbol x\in \R^n）。根据是否有其他的特殊约束，我们可以将连续优化再分为 无约束优化 (Unconstrained Optimization) 和 约束优化(Constrained Optimization) 两大类。

如，求 f(x)=x^2 的最小值显然是无约束优化问题；若限制 x\ge 1，则为约束优化问题。

无约束优化问题可以一般地表达为：

\min_{\boldsymbol x} f(\boldsymbol x)

约束优化问题表达为：[^3]

\begin{aligned} \min_{\boldsymbol x} \quad & f(\boldsymbol x)\\ \textrm{subject to} \quad & g_i(\boldsymbol x)=c_i & \textrm{for }i=1,\dots,n \quad & \textrm{Equality constraints} \\ \quad & h_j(\boldsymbol x)\ge c_j & \textrm{for }j=1,\dots,m \quad & \textrm{Inequality constraints} \end{aligned}

值得注意的是，线性规划问题是一类典型的约束优化问题。

[^3]: Constrained optimization - Wikipedia

梯度

在向量微积分中，梯度可简单解释为多元导数（微分）。对于关于 n 个变量 x_1,\dots,x_n 的函数 f : \R^n\to\R,\boldsymbol x\to f(\boldsymbol x),\boldsymbol x \in \R^n，我们分别对 f 求出关于 x_1,\dots,x_n 的偏导

\begin{aligned} \frac{\partial f}{\partial x_1}=&\lim_{h\to 0} \frac{f(x_1+h,x_2,\dots,x_n)-f(\boldsymbol x)}{h} \\ \vdots \\ \frac{\partial f}{\partial x_n}=&\lim_{h\to 0} \frac{f(x_1,\dots,x_{n-1},x_n+h)-f(\boldsymbol x)}{h} \\ \end{aligned}

并将它们合并到一个行向量中：

\nabla_x f=\mathrm{grad}f=\frac{\mathrm{d}f}{\mathrm{d}\boldsymbol x}=\begin{bmatrix} \frac{\partial f}{\partial x_1} & \frac{\partial f}{\partial x_2} & \cdots & \frac{\partial f}{\partial x_n} \end{bmatrix} \in \R^{1\times n}

这个行向量被称为 f 的 梯度 (Gradient) 。[^4]

[^4]: Mathematics for Machine Learning Draft 2024-01-15, Page 146

实际上，我们刚刚定义的梯度只是一种特殊情形。对于向量值域的函数 \boldsymbol f : \R^n\to\R^m 也有雅可比（Jacobian）矩阵 \boldsymbol J\in R^{m\times n}，可视为梯度的一般定义。本文不涉及向量值域函数，因此略过相关内容。感兴趣的读者可以自行阅读相关资料。[^5] [^6]

[^5]: 雅可比矩阵 - 维基百科，自由的百科全书 [^6]: Mathematics for Machine Learning Draft 2024-01-15, Page 150

梯度下降

梯度下降法的目标是解决一个多元函数 f : \R^n \to \R 的无约束优化问题：

\min_{\boldsymbol x} f(\boldsymbol x)

观察梯度的定义，不难发现函数上某点的梯度所指方向就是函数值增加最快的方向。因此，梯度的反方向即为函数值下降最快的方向。

如果我们从任意一个点开始，迭代地向梯度的反方向移动，那么最终一定能收敛到一个局部极小值。这一过程可以理解为一个小球在函数表面上自然滚动，直至停止。

具体如何实现呢？

我们选取一个很小的步长 \gamma，反复迭代 \boldsymbol x \larr \gamma (\nabla f)(\boldsymbol x)^\top。注意由于 \boldsymbol x 是列向量而梯度 (\nabla f)(\boldsymbol x) 是行向量，所以需要转置操作 ^\top。

对于合适的 \gamma，必然有 f(\boldsymbol x_0)\ge f(\boldsymbol x_1)\ge \dots 直至收敛。[^7]

在机器学习中，步长 \gamma 常常被称为 学习率 (Learning rate) 。[^8]

[^7]: Mathematics for Machine Learning Draft 2024-01-15, Page 228

[^8]: Learning rate - Wikipedia

我们来看一个具体的例子。考虑求 f(x)=x^4+7x^3+5x^2-17x+3 的最小值。其函数图像见下图（箭头表示梯度反方向）：[^9]

由于 f 是单变量函数，我们直接求导得

\frac{\mathrm{d}f}{\mathrm{d}x}=4x^3+21x^2+10x-17

根据基本常识，令导数等于 0 可解出三个顶点对应的 x：

我们希望求得全局最小值，即让 x 收敛到 -4.5 附近。若从 x=-3 处开始，沿着梯度反方向很快就能达到目标点；但如果从 x=-1 处开始，我们会在 x=1.6 的局部最小值附近不断游荡，无法达到目标。因此，x 的初始值的选择显得尤为重要。

如果目标函数是 凸函数 (Convex Function) ，显然从任意点开始迭代都能收敛到全局最小值。[^10]

[^9]: Mathematics for Machine Learning Draft 2024-01-15, Page 227 [^10]: 在 OI 中，凸函数的优化问题常常都能用三分法解决。介绍梯度下降，旨在提供一种常数更优、效率更高的做法。同时，梯度下降法可以并行优化多个变量，而嵌套三分会导致 \log 叠加。后面的例题中，我们将证实这一点。

How OI?

说了这么多，让我们先小试牛刀——使用梯度下降法通过弱化版三分板子。

例题 1：P3382 三分

题目大意：已知多项式函数 f(x) 在 [l,r] 上单峰。求

\mathop{\arg \max}_{x\in[l,r]} \ f(x)

即：取到 [l,r] 内最大值的 x。

考虑到多项式函数很容易求导，于是我们使用梯度下降法。这题是求 \max，我们可以改用“梯度上升”，即沿着梯度所指正方向移动。也可以令 g(x)=-f(x)，对 g 应用梯度下降。下面采用前一种做法。

由于本题限制了 x\in[l,r]，我们使用一个简单的改进措施：若一次迭代后 x 超出了 [l,r] 的区间，则撤销这次迭代，并减小步长。

我们先使用固定的初始步长和迭代次数，同时令 x_0=(l+r)/2。这里我使用了迭代次数 T=1000、步长 \gamma=0.01：

#include <cstdio>
using namespace std;

double a[15];

int main() {
    int n;
    double l, r;
    scanf("%d%lf%lf", &n, &l, &r);
    for(int i=n; i>=0; i--)
        scanf("%lf", a + i);
    for(int i=1; i<=n; i++)
        a[i - 1] = a[i] * i;
    auto diff = [&](double x) -> double {
        double res = 0;
        for(int i=n-1; i>=0; i--)
            (res *= x) += a[i];
        return res;
    };
    double x = (l + r) * 0.5;
    auto SGD = [&](int T, double lr) -> void {
        while(T--) {
            double grad = diff(x);
            double nx = x + grad * lr;
            if(nx < l || nx > r) lr *= 0.999;
            else x = nx;
        }
    };
    SGD(1000, 0.01);
    printf("%.10lf\n", x);
    return 0;
}

提交上去，喜提 33pts 的好成绩。

把迭代次数加到 10000？不好意思，44pts。

实际上，上面这份代码只修改 T 和 lr 的初始值是可以通过的。你可以自己玩玩 (?)

说到这里，你可能觉得梯度下降是一种玄学算法。没错，确实是。 不过，它的玄学程度远不如模拟退火，因为我们可以用一种简单直观的启发式方法来自适应步长：

我们具体化上面的步骤：减小步长即乘以 0.5，增大步长即乘以 1.1。同时，我们不再使用固定的迭代次数，改为当步长小于 \epsilon=10^{-7} 时停止。这样一来，狗都知道这算法不玄学了。下面是 100pts 的示例代码，建议先自己实现一遍再看。

#include <cstdio>
#define EPS 1e-7
using namespace std;

int n;
double a[15], b[15];

inline double eval(double x) {
    double res = 0;
    for(int i=0; i<=n; i++)
        (res *= x) += a[i];
    return res;
}

inline double diff(double x) {
    double res = 0;
    for(int i=0; i<n; i++)
        (res *= x) += b[i];
    return res;
}

int main() {
    double l, r;
    scanf("%d%lf%lf", &n, &l, &r);
    for(int i=0; i<=n; i++)
        scanf("%lf", a + i);
    for(int i=0; i<n; i++)
        b[i] = a[i] * (n - i);
    double x = (l + r) * 0.5;
    double lr = 1.0, grad = diff(x), val = eval(x);
    while(lr > EPS) {
        double nx = x + grad * lr;
        if(nx < l || nx > r) lr *= 0.8;
        else {
            double nval = eval(nx);
            if(nval > val)
                x = nx, grad = diff(nx), val = nval, lr *= 1.1;
            else lr *= 0.5;
        }
    }
    printf("%.10lf\n", x);
    return 0;
}

可喜的是，这个算法对超参数不是很敏感。如果你把 (1.1,0.5) 改成 (1.2,0.3) 甚至一些更激进的策略，经测试都能通过此题。初始步长也可以修改。在后面例题的实现中，我们将全部使用 (1.1,0.5,1.0) 这组超参数，\epsilon 则根据题目精度要求动态调整。

例题 2：P1337 平衡点

题意：求 n 个点的带权类费马点。

众所周知，这是一道模拟退火的板题。而我们将用梯度下降法来乱搞秒杀这道题，并通过运行时间证明梯度下降法在连续多元凸优化问题上的无比优越。（滑稽）

回归正题，我们要优化的目标函数为

f(x,y)=\sum_{i=1}^n w_i\sqrt{(x-x_i)^2+(y-y_i)^2}

求出 f 关于 x 和 y 的偏导分别为

\frac{\partial f}{\partial x}=\sum_{i=1}^n\frac{w_i(x-x_i)}{\sqrt{(x-x_i)^2+(y-y_i)^2}} \\ \frac{\partial f}{\partial y}=\sum_{i=1}^n\frac{w_i(y-y_i)}{\sqrt{(x-x_i)^2+(y-y_i)^2}}

于是，根据前面的自适应步长的思路可以写出代码：

// The code is left as an exercise for the reader.

这里其实有个很有意思的点。梯度的反方向实际上就是「合力」的方向。也就是说，我们本质上是在模拟物体的受力移动。这再次揭示了梯度下降法的本质：模拟轻质小球在曲面上滚动的过程。受到物理学的启发，科学家提出了 Momentum 优化算法，其基本思想即为模拟一个有质量的小球，在原版梯度下降的基础上添加「惯性」。限于篇幅，我们这里不介绍 Momentum，感兴趣的读者可以自行阅读 [^11][^12]。

[^11]: 12.6. Momentum — Dive into Deep Learning 1.0.3 documentation [^12]: Why Momentum Really Works

例题 3：P1883【模板】三分 | 函数

第一题太弱了，我们来试试真正的三分板子。

题意：有 n 个形如 ax^2+bx+c 的二次函数 f_1,f_2,\dots,f_n。设 F(x)=\max\{f_1(x),f_2(x),\dots,f_n(x)\}，求 F(x) 在 [0,1000] 上的最小值。

这道题其实跟第一题本质相同，不过计算 \max 的微分需要动点脑筋。

$$ \max'(f(x),g(x))=\begin{cases} f'(x) & f(x)>g(x) \\ g'(x) & f(x)<g(x) \end{cases} $$ 问题来了：当 $f(x)=g(x)$ 时，$\max$ 操作实际上是不可导的，因为此时从 $x$ 左边和右边逼近会得到不同的结果。 这时怎么办呢？其实我们任取 $f'(x)$ 和 $g'(x)$ 之一即可。在实际计算中，$f(x)$ 严格等于 $g(x)$ 的概率极小，以至于我们可以忽略这种情况。 > 上述处理方法参考了深度学习中常用的线性整流函数 $\mathrm{ReLU}$。其定义为 > $$ > \mathrm{ReLU}(x)=\max(0,x) \\ > \frac {\mathrm{d} \mathrm{ReLU}(x)}{\mathrm{d}x}=\begin{cases} > 0 & x<0 \\ > 1 & x>0 > \end{cases} > $$ > $\mathrm{ReLU}$ 函数在 $x=0$ 处不可微。实际操作中，我们将此时的梯度规定为 $0$ 或 $1$ 以保证训练的正常进行。 于是，我们得到： $$ F'(x)=f_k'(x) $$ 其中 $$ k=\mathop{\arg \max}_{i=1}^n \ f_i(x) $$ 即可。根据精度要求，取 $\epsilon=10^{-12}$。 ```cpp #include <cstdio> #define maxn 10005 #define EPS 1e-12 using namespace std; int n, a[maxn], b[maxn], c[maxn]; inline double eval(int i, double x) { return a[i] * x * x + b[i] * x + c[i]; } inline double diff(int i, double x) { return 2 * a[i] * x + b[i]; } inline int argmax(double x) { int best = 0; double best_val = eval(0, x); for(int i=1; i<n; i++) { double cur = eval(i, x); if(cur > best_val) best_val = cur, best = i; } return best; } const double l = 0, r = 1000; int main() { int T; scanf("%d", &T); while(T--) { scanf("%d", &n); for(int i=0; i<n; i++) scanf("%d%d%d", a + i, b + i, c + i); double x = (l + r) * 0.5; double lr = 1.0; int am = argmax(x); while(lr > EPS) { double nx = x - lr * diff(am, x); if(nx < l || nx > r) lr *= 0.8; else { int nm = argmax(nx); if(eval(nm, nx) < eval(am, x)) am = nm, x = nx, lr *= 1.1; else lr *= 0.5; } } printf("%.4lf\n", eval(am, x)); } return 0; } ``` 这份代码的超参数没有经过调优，基本上都是随手设的，也能跑到 $136\mathrm{ms}$，最优解第一页。还有很大优化空间！ ### 例题 4：[P4035 球形空间产生器](https://www.luogu.com.cn/problem/P4035) 题意：给定一个 $n$ 维球面上 $n+1$ 个点的坐标，请找出球心。$1\le n\le 10$。 我们先令这 $n+1$ 个坐标分别为 $\boldsymbol a_i=(a_{i,1},a_{i,2},\dots,a_{i,n})$（$0\le i\le n$）。 转换成数学语言，我们要找到下列关于 $\boldsymbol x\in \R^n$ 的方程组的一个近似解： $$ ||\boldsymbol x-\boldsymbol a_0||=||\boldsymbol x-\boldsymbol a_1||=\cdots=||\boldsymbol x-\boldsymbol a_n|| $$ 显然可以用高斯消元解决。不过这里我们介绍使用梯度下降的方法。 这里使用梯度下降的难点就在于目标函数的定义。如果我们定义目标函数 $$ f(x)=\begin{cases} 0 & \text{above equation holds} \\ 1 & \text{otherwise} \end{cases} $$ 则目标函数不连续，无法使用梯度下降优化。 我们换一种角度考虑这个问题。我们想让序列 $A=(||\boldsymbol x-\boldsymbol a_0||,\dots,||\boldsymbol x-\boldsymbol a_n||)$ 中的 $n+1$ 个元素**尽可能接近**。如何形式化表示？很容易想到**方差**。方差越小，说明我们的近似解 $\boldsymbol x$ 越接近精确解。当方差逼近 $0$ 时，我们可以近似认为这 $n+1$ 个实数是相等的。这样定义的目标函数在机器学习中称为**损失函数**。 最终我们的损失函数定义为**当前点到每个点距离的方差**，即 $$ f(\boldsymbol x)=\frac 1n\sum_{i=0}^n (||\boldsymbol x-\boldsymbol a_i||-b)^2 $$ 其中 $b$ 表示序列 $A$ 的算术平均数，即 $$ b=\frac 1n\sum_{i=0}^n ||\boldsymbol x-\boldsymbol a_i|| $$ 这个函数很难直接求导，一种解决方案是使用 **数值微分 (Numerical differentiation)** ，即令截距 $h$ 取一个非常小的正数（如 $10^{-9}$），计算 $$ \frac{\partial f}{\partial x} \approx \frac {f(x+h)-f(x)} {h} $$ 示例代码中使用**对称差分**： $$ \frac{\partial f}{\partial x} \approx \frac {f(x+h)-f(x-h)} {2h} $$ ```cpp #include <cstdio> #include <vector> #define EPS 1e-7 using namespace std; vector<double> p[11]; inline double sqr(double x) { return x * x; } inline double dist(const vector<double>& p1, const vector<double>& p2) { double sum = 0.0; for(int i=0; i<p1.size(); i++) sum += sqr(p1[i] - p2[i]); return __builtin_sqrt(sum); } inline double mean(const vector<double>& v) { double sum = 0.0; for(double x: v) sum += x; return sum / v.size(); } inline double variance(const vector<double>& v) { double m = mean(v); double sum = 0.0; for(double x: v) sum += sqr(x - m); return sum / v.size(); } inline double loss(const vector<double>& point) { int n = point.size(); vector<double> d(n + 1); for(int i=0; i<=n; i++) d[i] = dist(point, p[i]); return variance(d); } inline vector<double> numerical_gradient(vector<double>& point) { const double h = EPS; int n = point.size(); vector<double> grad(n); for(int i=0; i<n; i++) { double old = point[i]; point[i] += h; double loss_p = loss(point); point[i] = old - h; double loss_m = loss(point); point[i] = old; grad[i] = (loss_p - loss_m) / (2 * h); } return grad; } int main() { int n; scanf("%d", &n); for(int i=0; i<=n; i++) { p[i].resize(n); for(double &x: p[i]) scanf("%lf", &x); } vector<double> x(n, 0.0); double lr = 1.0; while(lr > EPS) { vector<double> grad = numerical_gradient(x); vector<double> x_new = x; for(int i=0; i<n; i++) x_new[i] -= lr * grad[i]; if(loss(x_new) < loss(x)) { x.swap(x_new); lr *= 1.1; } else { lr *= 0.5; } } for(double coord: x) printf("%.3lf ", coord); return 0; } ``` 另外也有深度学习中常用的基于**计算图**[^13]的方法。它对目标函数的完整运算过程构建了一个 DAG，并应用**链式求导法则**反向传播梯度。有兴趣的读者可以自行尝试。 [^13]: 原书：[深度学习入门：基于Python的理论与实现](https://www.ituring.com.cn/book/1921) 搬运博文：[机器学习入门（10）— 浅显易懂的计算图、链式法则讲解](https://blog.csdn.net/wohu1104/article/details/106911071) ## 算法对比 & 局限性 下面的表格展示了梯度下降、三分、模拟退火三种算法各自的优势： | | 梯度下降 | 三分 | 模拟退火 | | :--------: | :-------------------: | :--------------------------: | :-------------------: | | 连续优化 | ✅ | ✅ | ❔ [^14] | | 离散优化 | ❌ | ✅ | ✅ | | 非凸优化 | ❔ [^15] | ❌ | ✅ | | 时间复杂度 | 玄学，<三分 | $\mathcal O(\log^k V)$ [^16] | 玄学，>三分 | [^14]: 模拟退火的精度难以控制。 [^15]: 梯度下降算法收敛到局部还是全局最小值，取决于超参数（初始值、步长）的设置。也可以随机多个点分别跑梯度下降来提升成功率。 [^16]: $k$ 表示自变量个数，$V$ 表示值域规模。 ## 总结 好消息！好消息！全新的乱搞算法——梯度下降出炉啦！！ 不要 TLE，不要 $1.5\mathrm s$，只要 $3\mathrm{ms}$！ 微分大法好，梯度保平安。 左脚踹翻臭三分，右拳暴击 S(imulate) A(nnealing)。 连续优化都除尽，GD 一统靖江山。 ## References & Notes