[最优化技术] 3-9 共轭梯度法

3-9 共轭梯度法

共轭梯度法

简述

\(\qquad\)上一篇介绍了共轭方向法。它用平行搜索构造共轭方向：从两个点出发沿同一方向各搜索一次，把两个一维极小点连起来，就得到与搜索方向共轭的新方向。这个方法只需要函数值，但每构造一个新方向都要作两次一维搜索；要用它系统地得到 \(n\) 个共轭方向，还得配合方向替换。

\(\qquad\)上一篇末尾已经埋下伏笔：二次函数的梯度差满足 \(\nabla f(X^{(k+1)}) - \nabla f(X^{(k)}) = H(X^{(k+1)} - X^{(k)})\)，所以共轭条件里的矩阵 \(H\) 可以用相邻两点的梯度差代替。只要利用迭代点处的梯度，就能直接递推出共轭方向，既不必计算 \(H\)，也不必为了构造一个方向去作两次平行搜索。这就是共轭梯度法（Conjugate Gradient Method）。

\(\qquad\)共轭梯度法是共轭方向法中的一种，因为它的每个共轭向量都依赖于迭代点处的负梯度构造出来。它的第一个搜索方向取初始点的负梯度方向，这就是最速下降法；以后各步的搜索方向是把负梯度偏转一个角度，相当于对负梯度作了修正。所以共轭梯度法实质上是对最速下降法的改进，又称旋转梯度法。

共轭方向与梯度的关系

考虑正定二次函数

\[f(X) = \frac{1}{2}X^T H X + B^T X + C \]

在点 \(X^{(k)}\)、\(X^{(k+1)}\) 处的梯度分别为

\[\nabla f(X^{(k)}) = HX^{(k)} + B, \qquad \nabla f(X^{(k+1)}) = HX^{(k+1)} + B \]

两式相减，得

\[\nabla f(X^{(k+1)}) - \nabla f(X^{(k)}) = H\left(X^{(k+1)} - X^{(k)}\right) \]

若 \(X^{(k+1)}\) 是从 \(X^{(k)}\) 沿 \(S^{(k)}\) 作一维搜索得到的，则 \(X^{(k+1)} - X^{(k)} = \alpha^{(k)}S^{(k)}\)，上式写成

\[\nabla f(X^{(k+1)}) - \nabla f(X^{(k)}) = \alpha^{(k)} H S^{(k)} \]

用前面已经得到的某个共轭方向 \(S^{(j)}\ (j < k)\) 左乘上式：

\[[S^{(j)}]^T\left[\nabla f(X^{(k+1)}) - \nabla f(X^{(k)})\right] = \alpha^{(k)}[S^{(j)}]^T H S^{(k)} \]

而 \(S^{(j)}\) 与 \(S^{(k)}\) 共轭，右端为零，于是

\[[S^{(j)}]^T\left[\nabla f(X^{(k+1)}) - \nabla f(X^{(k)})\right] = 0 \]

\(\qquad\)这就是共轭方向与梯度之间的关系：相邻两点的梯度差，与所有已用的共轭方向都正交。反过来，只要梯度差与 \(S^{(j)}\) 正交，由于 \(\alpha^{(k)} \ne 0\)，\(S^{(j)}\) 与 \(S^{(k)}\) 必然关于 \(H\) 共轭，两者是等价的。另一方面，\(X^{(k+1)}\) 是沿 \(S^{(k)}\) 的一维极小点，所以 \([S^{(k)}]^T\nabla f(X^{(k+1)}) = 0\)。两条合起来，共轭条件可以完全用梯度表达，\(H\) 不必出现——这正是共轭梯度法能够只靠梯度构造方向的原因。

\(\qquad\)图 4-10 画出了这个关系：\(g_k\)、\(g_{k+1}\) 分别是 \(X^k\)、\(X^{k+1}\) 处的梯度，\(g_{k+1} - g_k\) 垂直于已有的共轭方向 \(d^j\)，而 \(g_{k+1}\) 又垂直于本次搜索方向 \(d^k\)。

符号说明：图中 \(x^k\)、\(x^{k+1}\) 即正文的 \(X^{(k)}\)、\(X^{(k+1)}\)；\(g_k\)、\(g_{k+1}\) 即 \(\nabla f(X^{(k)})\)、\(\nabla f(X^{(k+1)})\)；\(d^j\) 是已得到的共轭方向 \(S^{(j)}\)，\(d^k\) 是本次的搜索方向 \(S^{(k)}\)。

递推公式

既然共轭条件可以用梯度差表达，就可以在负梯度方向上作一次修正，构造出与上一个方向共轭的新方向。取

\[S^{(k+1)} = -\nabla f(X^{(k+1)}) + \beta^{(k)} S^{(k)} \]

即在新点的负梯度上，沿上一步的搜索方向偏转一个由 \(\beta^{(k)}\) 控制的份额。系数 \(\beta^{(k)}\) 由共轭条件 \([S^{(k+1)}]^T H S^{(k)} = 0\) 确定。

先看第一步。记 \(g_0 = \nabla f(X^{(0)})\)、\(g_1 = \nabla f(X^{(1)})\)，此时 \(S^{(0)} = -g_0\)。由

\[g_1 - g_0 = \alpha^{(0)} H S^{(0)} \]

把 \(S^{(1)} = -g_1 + \beta^{(0)}S^{(0)}\) 代入共轭条件，得

\[[S^{(1)}]^T H S^{(0)} = \frac{1}{\alpha^{(0)}}\left[-g_1 + \beta^{(0)}S^{(0)}\right]^T\left(g_1 - g_0\right) = 0 \]

展开，注意 \(S^{(0)} = -g_0\)，且 \(g_1\) 与 \(S^{(0)}\) 正交（一维搜索），因而 \(g_1\) 与 \(g_0\) 也正交，于是

\[-\|g_1\|^2 + \beta^{(0)}\|g_0\|^2 = 0 \]

即

\[\beta^{(0)} = \frac{\|g_1\|^2}{\|g_0\|^2} = \frac{\|\nabla f(X^{(1)})\|^2}{\|\nabla f(X^{(0)})\|^2} \]

对第 \(k\) 步重复同样的过程，得到一般的递推式

\[\beta^{(k)} = \frac{\|\nabla f(X^{(k+1)})\|^2}{\|\nabla f(X^{(k)})\|^2} \]

\[S^{(k+1)} = -\nabla f(X^{(k+1)}) + \beta^{(k)}S^{(k)} \]

\(\qquad\)这个系数只用相邻两点的梯度模，不涉及 \(H\)，也不涉及任何矩阵运算。这就是共轭梯度法与一般共轭方向法的区别：方向仍然关于 \(H\) 共轭，但构造过程只用到梯度。

\(\qquad\)注意：有的资料把第二项写成 \(-\beta^{(k)}\nabla f(X^{(k)})\)。由于第一步 \(S^{(0)} = -\nabla f(X^{(0)})\)，两种写法在第一步完全一致；但递推时被继承下来的是上一步的搜索方向，所以一般写成 \(S^{(k)}\)。

具体步骤

对于给定初始点 \(X^{(0)}\) 和收敛精度 \(\varepsilon\)：

① 令 \(k = 0\)；

② 取初始点的负梯度方向为第一个搜索方向

\[S^{(0)} = -\nabla f(X^{(0)}) \]

③ 沿 \(S^{(k)}\) 作一维搜索，求最优步长 \(\alpha^{(k)}\)，使

\[f(X^{(k)} + \alpha^{(k)}S^{(k)}) = \min_{\alpha} f(X^{(k)} + \alpha S^{(k)}) \]

并令

\[X^{(k+1)} = X^{(k)} + \alpha^{(k)}S^{(k)} \]

④ 若 \(\|\nabla f(X^{(k+1)})\| < \varepsilon\)，则停止，输出 \(X^* = X^{(k+1)}\)；否则转⑤；

⑤ 若 \(k = n\)，则令 \(X^{(0)} = X^{(k+1)}\)，转②重新取负梯度方向（重置）；否则转⑥；

⑥ 计算

\[\beta^{(k)} = \frac{\|\nabla f(X^{(k+1)})\|^2}{\|\nabla f(X^{(k)})\|^2}, \qquad S^{(k+1)} = -\nabla f(X^{(k+1)}) + \beta^{(k)}S^{(k)} \]

令 \(k = k + 1\)，转③。

符号说明：

• \(X^{(k)}\) 是第 \(k\) 次迭代点，也是第 \(k-1\) 次搜索的终点；
• \(S^{(k)}\) 是第 \(k\) 个搜索方向，\(n\) 个方向关于 \(H\) 两两共轭；
• \(\alpha^{(k)}\) 是沿 \(S^{(k)}\) 的最优步长，\(\beta^{(k)}\) 是构造下一个方向的共轭系数。

\(\qquad\)第⑤步对二次函数来说通常不会走到：\(n\) 个共轭方向用完，\(n\) 元二次函数已经到达极小点。它的作用在非二次函数上——此时 \(n\) 次搜索后共轭性已经被函数的非二次项破坏，继续沿用旧方向没有意义，不如重新从负梯度出发，重新积累一组共轭方向。

程序框图

符号说明：图中 \(x\)、\(d^k\)、\(g_k\)、\(\alpha_k\)、\(\beta_k\) 分别对应正文的 \(X\)、\(S^{(k)}\)、\(\nabla f(X^{(k)})\)、\(\alpha^{(k)}\)、\(\beta^{(k)}\)；判据为 \(\|\nabla f(X^{(k+1)})\| < \varepsilon\)；“\(k = n\)”框即第⑤步的重置，重置后回到取 \(S^{(0)} = -\nabla f(X^{(0)})\) 的框。

实例计算

仍用前面几篇的标准函数，完整地算两步。

题目

\[f(X) = x_1^2 + 2x_2^2 - 4x_1 - 2x_1x_2,\quad X^{(0)} = [1,1]^T,\quad \varepsilon = 0.01 \]

它的梯度和 Hessian矩阵为

\[\nabla f(X) = \begin{bmatrix} 2x_1 - 2x_2 - 4 \\ 4x_2 - 2x_1 \end{bmatrix},\qquad H = \begin{bmatrix} 2 & -2 \\ -2 & 4 \end{bmatrix},\qquad B = \begin{bmatrix} -4 \\ 0 \end{bmatrix} \]

理论最优解为 \(X^* = [4,2]^T\)，\(f(X^*) = -8\)。

解答

1. 第一次搜索

\[\nabla f(X^{(0)}) = \begin{bmatrix} 2 - 2 - 4 \\ 4 - 2 \end{bmatrix} = \begin{bmatrix} -4 \\ 2 \end{bmatrix},\qquad S^{(0)} = -\nabla f(X^{(0)}) = \begin{bmatrix} 4 \\ -2 \end{bmatrix} \]

沿 \(S^{(0)}\) 搜索，此时 \(X = [1+4\alpha,\ 1-2\alpha]^T\)，

\[\begin{aligned} \varphi(\alpha) &= (1+4\alpha)^2 + 2(1-2\alpha)^2 - 4(1+4\alpha) - 2(1+4\alpha)(1-2\alpha) \\ &= 40\alpha^2 - 20\alpha - 3 \end{aligned} \]

令 \(\varphi'(\alpha) = 80\alpha - 20 = 0\)，得 \(\alpha^{(0)} = 0.25\)，所以

\[X^{(1)} = \begin{bmatrix} 1 \\ 1 \end{bmatrix} + 0.25\begin{bmatrix} 4 \\ -2 \end{bmatrix} = \begin{bmatrix} 2 \\ 0.5 \end{bmatrix},\qquad f(X^{(1)}) = -5.5 \]

该点的梯度为

\[\nabla f(X^{(1)}) = \begin{bmatrix} 4 - 1 - 4 \\ 2 - 4 \end{bmatrix} = \begin{bmatrix} -1 \\ -2 \end{bmatrix},\qquad \|\nabla f(X^{(1)})\| = \sqrt{(-1)^2 + (-2)^2} = \sqrt{5} \approx 2.236 > 0.01 \]

未达精度，继续。

2. 构造共轭方向

\[\beta^{(0)} = \frac{\|\nabla f(X^{(1)})\|^2}{\|\nabla f(X^{(0)})\|^2} = \frac{(-1)^2 + (-2)^2}{(-4)^2 + 2^2} = \frac{5}{20} = 0.25 \]

\[S^{(1)} = -\nabla f(X^{(1)}) + \beta^{(0)}S^{(0)} = \begin{bmatrix} 1 \\ 2 \end{bmatrix} + 0.25\begin{bmatrix} 4 \\ -2 \end{bmatrix} = \begin{bmatrix} 2 \\ 1.5 \end{bmatrix} \]

验算共轭性：

\[H S^{(1)} = \begin{bmatrix} 2 & -2 \\ -2 & 4 \end{bmatrix}\begin{bmatrix} 2 \\ 1.5 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \end{bmatrix},\qquad [S^{(0)}]^T H S^{(1)} = [4,-2]\begin{bmatrix} 1 \\ 2 \end{bmatrix} = 0 \]

确实共轭。

3. 第二次搜索

沿 \(S^{(1)}\) 搜索，此时 \(X = [2+2\alpha,\ 0.5+1.5\alpha]^T\)，

\[\begin{aligned} \varphi(\alpha) &= (2+2\alpha)^2 + 2(0.5+1.5\alpha)^2 - 4(2+2\alpha) - 2(2+2\alpha)(0.5+1.5\alpha) \\ &= 2.5\alpha^2 - 5\alpha - 5.5 \end{aligned} \]

令 \(\varphi'(\alpha) = 5\alpha - 5 = 0\)，得 \(\alpha^{(1)} = 1\)，所以

\[X^{(2)} = \begin{bmatrix} 2 \\ 0.5 \end{bmatrix} + 1\cdot\begin{bmatrix} 2 \\ 1.5 \end{bmatrix} = \begin{bmatrix} 4 \\ 2 \end{bmatrix} \]

此时

\[\nabla f(X^{(2)}) = \begin{bmatrix} 8 - 4 - 4 \\ 8 - 8 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix},\qquad f(X^{(2)}) = -8 \]

满足 \(\|\nabla f(X^{(2)})\| = 0 < \varepsilon\)，停止，\(X^* = [4,2]^T\)。

把两次搜索列成表：

次序	起点	\(\nabla f(X^{(k)})\)	\(\beta^{(k-1)}\)	\(S^{(k)}\)	\(\alpha^{(k)}\)	终点	\(f\)
1	\((1,1)\)	\([-4,2]^T\)	—	\([4,-2]^T\)	\(0.25\)	\((2,0.5)\)	\(-5.5\)
2	\((2,0.5)\)	\([-1,-2]^T\)	\(0.25\)	\([2,1.5]^T\)	\(1\)	\((4,2)\)	\(-8\)

\(\qquad\)两次搜索就到达精确极小点，正好是 \(n=2\) 步。第一步走的是最速下降方向，如果继续沿负梯度走，就会像最速下降法那样出现锯齿；共轭梯度法在第二步把负梯度 \([1,2]^T\) 偏转成 \([2,1.5]^T\)，一步到位。

与前两篇的方法对比：

方法	是否用梯度	一维搜索次数	结果	\(f\)
坐标轮换法（\(\varepsilon = 0.01\)）	否	\(8 \times 2 = 16\)	\((3.992188, 1.996094)\)	\(-7.999969\)
共轭方向法	否	\(3\)	\((4,2)\)	\(-8\)
共轭梯度法	是	\(2\)	\((4,2)\)	\(-8\)

代码示例

\(\qquad\)程序沿用进退法和黄金分割法作一维搜索。与上一篇相比，这里不需要第二个初始点，也不需要平行搜索；每次迭代只做三件事：沿当前方向一维搜索、算新点的梯度、用梯度模的比值构造下一个方向。方向用满 \(n\) 个以后，按第⑤步重置为负梯度方向。

核心的递推部分如下：

C++代码示例（共轭方向递推）：

        // 3. 构造新的共轭方向
        if ((k + 1) % n == 0) {
            // 已用满 n 个方向，重置为负梯度方向
            for (int i = 0; i < n; ++i) s[i] = -g_new[i];
            cout << "  k = n，重置 S = -grad f" << endl;
        } else {
            double g2_old = 0.0;
            for (double v : g) g2_old += v * v;
            double beta = g2 / g2_old;
            for (int i = 0; i < n; ++i) s[i] = -g_new[i] + beta * s[i];
            cout << "  beta = " << beta << ", ";
            print_vec("S", s);
            cout << endl;
        }

完整程序如下，目标函数、梯度、一维搜索和输出都在其中。

C++代码实现（完整）：

#include <iostream>
#include <vector>
#include <cmath>
#include <functional>
#include <iomanip>
#include <stdexcept>
#include <utility>

using namespace std;

// 进退法求初始区间 (复用)
pair<double, double> bracket_minimum(function<double(double)> func, double a0, double h) {
    const int MAX_ITER = 1e4;
    if (h <= 0) throw invalid_argument("Step size h must be positive.");
    if (func(a0) <= func(a0 + h)) {
        double x0 = a0 + h, x1 = a0;
        h *= 2;
        double x2 = x1 - h;
        int iter = 0;
        while (func(x2) < func(x1)) {
            if (iter++ > MAX_ITER) throw runtime_error("Max iterations exceeded");
            x0 = x1; x1 = x2; h *= 2; x2 = x1 - h;
        }
        return make_pair(x2, x0);
    } else {
        double x0 = a0;
        double x1 = a0 + h;
        h *= 2;
        double x2 = x1 + h;
        int iter = 0;
        while (func(x2) < func(x1)) {
            if (iter++ > MAX_ITER) throw runtime_error("Max iterations exceeded");
            x0 = x1; x1 = x2; h *= 2; x2 = x1 + h;
        }
        return make_pair(x0, x2);
    }
}

// 黄金分割法 (复用)
double golden_section(function<double(double)> func, double l, double r, double eps) {
    double a = l, b = r, aa, bb;
    while (b - a > eps) {
        aa = b - 0.618 * (b - a);
        bb = a + 0.618 * (b - a);
        if (func(aa) < func(bb)) b = bb;
        else a = aa;
    }
    return (a + b) / 2;
}

// 从点 x 出发，沿方向 s 作一维搜索
vector<double> line_search(
    function<double(const vector<double>&)> func,
    const vector<double>& x,
    const vector<double>& s,
    double eps
) {
    auto phi = [&](double alpha) {
        vector<double> y(x.size());
        for (size_t i = 0; i < x.size(); ++i) y[i] = x[i] + alpha * s[i];
        return func(y);
    };

    pair<double, double> bracket;
    try { bracket = bracket_minimum(phi, 0.0, 0.1); }
    catch (...) { bracket = make_pair(0.0, 10.0); }

    double alpha_star = golden_section(phi, bracket.first, bracket.second, eps);

    vector<double> y(x.size());
    for (size_t i = 0; i < x.size(); ++i) y[i] = x[i] + alpha_star * s[i];
    return y;
}

// 输出一个向量
void print_vec(const char* name, const vector<double>& v) {
    cout << name << " = [";
    for (size_t i = 0; i < v.size(); ++i) cout << v[i] << (i == v.size() - 1 ? "" : ", ");
    cout << "]^T";
}

// 共轭梯度法：f 为目标函数，grad 为梯度函数，x0 为初始点
vector<double> conjugate_gradient(
    function<double(const vector<double>&)> f,
    function<vector<double>(const vector<double>&)> grad,
    const vector<double>& x0,
    double eps
) {
    int n = x0.size();
    vector<double> x = x0;
    vector<double> g = grad(x);

    // S^(0) = -grad f(X^(0))
    vector<double> s(n);
    for (int i = 0; i < n; ++i) s[i] = -g[i];

    for (int k = 0; k < 1000; ++k) {
        // 1. 沿 S^(k) 作一维搜索，得到 X^(k+1)
        x = line_search(f, x, s, eps);
        vector<double> g_new = grad(x);

        double g2 = 0.0;
        for (double v : g_new) g2 += v * v;
        double norm_g = sqrt(g2);

        cout << "第 " << k + 1 << " 次搜索: ";
        print_vec("X", x);
        cout << ", ||grad f|| = " << norm_g << endl;

        // 2. 收敛判断
        if (norm_g < eps) {
            cout << "  满足 ||grad f|| < " << eps << "，停止迭代" << endl;
            break;
        }

        // 3. 构造新的共轭方向
        if ((k + 1) % n == 0) {
            // 已用满 n 个方向，重置为负梯度方向
            for (int i = 0; i < n; ++i) s[i] = -g_new[i];
            cout << "  k = n，重置 S = -grad f" << endl;
        } else {
            double g2_old = 0.0;
            for (double v : g) g2_old += v * v;
            double beta = g2 / g2_old;
            for (int i = 0; i < n; ++i) s[i] = -g_new[i] + beta * s[i];
            cout << "  beta = " << beta << ", ";
            print_vec("S", s);
            cout << endl;
        }

        g = g_new;
    }
    return x;
}

int main() {
    cout << fixed << setprecision(6);
    cout << "共轭梯度法求解二元二次函数极小值" << endl;
    cout << "目标函数: f(X) = x1^2 + 2*x2^2 - 4*x1 - 2*x1*x2" << endl;
    cout << "理论最优解: X* = [4.0, 2.0]^T, f(X*) = -8.0" << endl << endl;

    auto f = [](const vector<double>& x) {
        return x[0] * x[0] + 2 * x[1] * x[1] - 4 * x[0] - 2 * x[0] * x[1];
    };
    auto grad = [](const vector<double>& x) {
        return vector<double>{2 * x[0] - 2 * x[1] - 4, 4 * x[1] - 2 * x[0]};
    };

    vector<double> x0 = {1.0, 1.0};
    double eps = 0.01;

    vector<double> result = conjugate_gradient(f, grad, x0, eps);

    cout << "\n最终结果:" << endl;
    print_vec("最优解 X*", result);
    cout << endl;
    cout << "目标函数值 f(X*) = " << f(result) << endl;

    return 0;
}

\(\qquad\)程序的实测过程如下：第一次搜索得到 \((2.001808, 0.499096)\)，\(\|\nabla f\| = 2.240126\)，算得 \(\beta^{(0)} = 0.250908\)，方向为 \([1.998208, 1.505417]^T\)；第二次搜索得到 \((3.990933, 1.997670)\)，此时 \(\|\nabla f\| = 0.016100\)，还没有达到 \(0.01\) 的精度，而两个方向已经用完，于是按第⑤步重置为负梯度方向，再搜一次，得到 \(X^* = (3.993963, 1.995688)\)，\(f = -7.999978\)。

\(\qquad\)手算两步就到 \((4,2)\)，是因为一维搜索取了精确的极小点；程序用黄金分割法，精度 \(0.01\)，所以第二步停在 \((3.990933, 1.997670)\)，与理论值相差约 \(0.009\)，最后的误差就来自这里。

方法分析

1. 只需要梯度：每步作一次一维搜索、算一次梯度，不形成也不求逆 Hessian矩阵。存储量就是几个 \(n\) 维向量，与最速下降法相同，比牛顿法省得多。

2. \(n\) 步收敛：对 \(n\) 元正定二次函数，只要一维搜索是精确的，\(n\) 个共轭方向用完就能到达极小点。这是共轭方向的性质在梯度构造下的体现。

3. 是对最速下降法的改进：第一步与最速下降法相同，以后每步把负梯度偏转一个角度，避免相邻方向正交造成的锯齿，所以又称旋转梯度法。方向一次比一次“有效”，不像坐标轮换法那样沿固定的坐标轴方向反复搜索。

4. 与共轭方向法的区别在构造方式：平行搜索每构造一个方向要作两次一维搜索，还要另选一个初始点；共轭梯度法用梯度递推，一次搜索就前进一个方向，代价是要能计算梯度。

5. 共轭性依赖一维搜索的精度：一维搜索不精确，方向的共轭性就只是近似的，\(n\) 步收敛也不再严格。非二次函数的情形更是如此，所以算法保留了第⑤步的重置：用满 \(n\) 个方向就重新从负梯度出发，重新积累共轭方向。

总结

\(\qquad\)共轭梯度法把上一篇的“构造共轭方向”换了一个做法：不再用两次平行搜索去凑，而是用相邻两点的梯度差把 \(H\) 消掉，直接在负梯度上作修正。方向仍然关于 \(H\) 共轭，二次函数仍然 \(n\) 步收敛，但每一步只需要一次一维搜索和一次梯度计算。

\(\qquad\)它的代价也很清楚：必须能算梯度。下一篇介绍鲍威尔法，它回到只用函数值的路子上，用平行搜索加方向替换来组织共轭方向，不需要梯度，也不需要 \(H\)，是直接搜索法里很有代表性的一种。

---

本学习笔记参考资料：

[1] 白清顺, 孙靖民, 梁迎春. 机械优化设计 第7版[M]. 北京: 机械工业出版社, 2024. ISBN: 978-7-111-75103-8 在线链接

[2] 武汉理工大学《最优化技术B》课程课件，授课教师：颜彬老师