[最优化技术] 3-5 牛顿法与阻尼牛顿法

3-5 牛顿法与阻尼牛顿法

牛顿法

简述

\(\qquad\)牛顿法（Newton's Method）是一种利用目标函数的二阶导数信息（Hessian矩阵）进行优化的方法。与最速下降法仅使用一阶梯度信息不同，牛顿法通过二次泰勒展开近似原函数，能够更准确地预测极值点的位置，因此具有更快的收敛速度（二阶收敛）。

牛顿法的基本思想是：在当前点处用二次函数近似原目标函数，然后直接求出该二次函数的极小值点作为下一次迭代点。

原理

对于多元函数 \(f(X)\)，在点 \(X^{(k)}\) 处进行二阶泰勒展开：

\[f(X) \approx f(X^{(k)}) + [\nabla f(X^{(k)})]^T (X - X^{(k)}) + \frac{1}{2}(X - X^{(k)})^T H(X^{(k)}) (X - X^{(k)}) \]

其中 \(H(X^{(k)})\) 为 Hessian矩阵（海森矩阵）：

\[H(f) = \left(\begin{array}{cccc} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_1 \partial x_n} \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} & \cdots & \frac{\partial^2 f}{\partial x_2 \partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2 f}{\partial x_n \partial x_1} & \frac{\partial^2 f}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_n^2} \end{array}\right) \]

注意：当函数的所有二阶混合偏导数连续时，有 \(\frac{\partial^2 f}{\partial x_i \partial x_j} = \frac{\partial^2 f}{\partial x_j \partial x_i}\)，即 H矩阵对称。对于多项式函数，Hessian矩阵必然对称。

对于二次函数求极值，令泰勒展开的一阶导数为零：

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

解得牛顿法迭代公式：

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

Hessian矩阵的性质

Hessian矩阵的正定性决定了极值点的性质：

  1. 正定（所有特征值 > 0）→ 局部极小点
  2. 负定（所有特征值 < 0）→ 局部极大点
  3. 不定（有正有负特征值）→ 鞍点
  4. 半定（有零特征值）→ 需要高阶检验

对于二元二次函数，Hessian矩阵为：

\[H = \left(\begin{array}{cc} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} \end{array}\right) \]

二阶矩阵求逆公式：对于矩阵 \(A = \left(\begin{array}{cc} a & b \\ c & d \end{array}\right)\)，当 \(ad - bc \neq 0\) 时：

\[A^{-1} = \frac{1}{ad-bc} \left(\begin{array}{cc} d & -b \\ -c & a \end{array}\right) \]

牛顿法分析

特点：

  1. 原函数为二次函数：一步求出极值点。
  2. 原函数为非二次函数：先通过二阶泰勒展开变为二次函数，因此求出的为极值点的近似值。

优点：

  1. 收敛速度最快（二阶收敛）
  2. 对于二次函数，一次迭代即可收敛

缺点：

  1. Hessian矩阵计算困难（需要计算所有二阶偏导数）
  2. Hessian矩阵求逆计算量大
  3. 当Hessian矩阵不正定时，可能不收敛或收敛到鞍点

具体步骤

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

  • ① 确定初始点 \(X^{(0)}\)，收敛精度 \(\varepsilon\)，令 \(k=0\)；
  • ② 计算梯度 \(\nabla f(X^{(k)})\) 和 Hessian矩阵 \(H(X^{(k)})\)；
  • ③ 判断Hessian矩阵是否正定。若不正定，需采用修正策略（详见下一篇文章《3-6 修正牛顿法》）；
  • ④ 构造搜索方向 \(d^{(k)} = -[H(X^{(k)})]^{-1} \nabla f(X^{(k)})\)；
  • ⑤ 更新迭代点：\(X^{(k+1)} = X^{(k)} + d^{(k)}\)，\(k=k+1\)；
  • ⑥ 重复②~⑤步，直到：\(\|\nabla f(X^{(k)})\| \le \varepsilon\)。

实例计算

这里我们举个例子，手把手带大家使用牛顿法，一步步来计算一个二元二次凸函数的最优解。

题目

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

解答

首先计算 \(f(X)\) 的梯度函数和Hessian矩阵：

梯度：

\[\nabla f(X) = \left(\begin{array}{c} 2x_1 - 2x_2 -4 \\ -2x_1 + 4x_2 \end{array}\right) \]

Hessian矩阵：

\[H(X) = \left(\begin{array}{cc} \frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} \end{array}\right) = \left(\begin{array}{cc} 2 & -2 \\ -2 & 4 \end{array}\right) \]

注意：对于二次函数，Hessian矩阵为常数矩阵，不随 \(X\) 变化。

代入 \(X^{(0)} = [1,1]^T\)，得：

\[\nabla f(X^{(0)}) = \left(\begin{array}{c} 2(1) - 2(1) -4 \\ -2(1) + 4(1) \end{array}\right) = \left(\begin{array}{c} -4 \\ 2 \end{array}\right) \]

计算梯度的模长：

\[\|\nabla f(X^{(0)})\| = \sqrt{(-4)^2 + 2^2} = \sqrt{20} = 2\sqrt{5} \approx 4.47 > \varepsilon \]

不满足收敛条件，迭代继续。

第 1 次迭代 (\(k=0\))：

  1. 计算Hessian矩阵的逆：

\[H = \left(\begin{array}{cc} 2 & -2 \\ -2 & 4 \end{array}\right) \]

行列式：\(\det(H) = 2 \times 4 - (-2) \times (-2) = 8 - 4 = 4\)

\[H^{-1} = \frac{1}{4} \left(\begin{array}{cc} 4 & 2 \\ 2 & 2 \end{array}\right) = \left(\begin{array}{cc} 1 & 0.5 \\ 0.5 & 0.5 \end{array}\right) \]

  1. 计算牛顿方向：

\[d^{(0)} = -H^{-1} \nabla f(X^{(0)}) = -\left(\begin{array}{cc} 1 & 0.5 \\ 0.5 & 0.5 \end{array}\right) \left(\begin{array}{c} -4 \\ 2 \end{array}\right) = -\left(\begin{array}{c} -4 + 1 \\ -2 + 1 \end{array}\right) = \left(\begin{array}{c} 3 \\ 1 \end{array}\right) \]

  1. 更新迭代点：

\[X^{(1)} = X^{(0)} + d^{(0)} = \left(\begin{array}{c} 1 \\ 1 \end{array}\right) + \left(\begin{array}{c} 3 \\ 1 \end{array}\right) = \left(\begin{array}{c} 4 \\ 2 \end{array}\right) \]

  1. 检验收敛条件：
    计算 \(X^{(1)}\) 处的梯度：

\[\nabla f(X^{(1)}) = \left(\begin{array}{c} 2(4) - 2(2) - 4 \\ -2(4) + 4(2) \end{array}\right) = \left(\begin{array}{c} 0 \\ 0 \end{array}\right) \]

梯度的模长：

\[\|\nabla f(X^{(1)})\| = 0 < \varepsilon \]

满足收敛条件，迭代结束！

最终得到最优解：\(X^* = X^{(1)} = [4, 2]^T\)，\(f(X^*) = 4^2 + 2 \times 2^2 - 4 \times 4 - 2 \times 4 \times 2 = 16 + 8 - 16 - 16 = -8\)

结论：当目标函数 \(f(X)\) 为二次函数时，二阶泰勒展开是精确的，而其中的Hessian矩阵 \(\nabla^2f(X^{(k)})\) 是一个常数矩阵。因此从任意初始点进行迭代，只需一步迭代即可找到目标函数的极小值点。

代码示例

牛顿法 C++ 代码示例： 
// 牛顿法 求多元函数最小值点（支持任意维度）
// func: 目标函数
// grad: 梯度函数
// hess: Hessian矩阵函数
// x0: 初始点（向量）
// eps: 收敛精度
vector<double> newton_method(
    function<double(const vector<double>&)> func,
    function<vector<double>(const vector<double>&)> grad,
    function<vector<vector<double>>(const vector<double>&)> hess,
    const vector<double>& x0,
    double eps,
    int max_iter = 1000
) {
    vector<double> x = x0;
    int n = x0.size();
    
    cout << "\n迭代过程:" << endl;
    cout << "步骤\t点 X\t\t\t函数值 f(X)\t\t梯度范数" << endl;
    cout << "-----------------------------------------------------------------" << endl;
    
    for (int iter = 0; iter < max_iter; ++iter) {
        // 1. 计算当前点的梯度
        vector<double> g = grad(x);
        
        // 2. 计算梯度范数 ||∇f(X)||
        double grad_norm = norm(g);
        
        // 3. 检查收敛条件
        if (grad_norm < eps) {
            cout << iter << "\t";
            for (size_t i = 0; i < x.size(); ++i) {
                printf("%.6lf", x[i]);
                cout << (i == x.size() - 1 ? "" : ", ");
            }
            printf("\t%.6lf\t\t%.6lf (收敛)\n", func(x), grad_norm);
            break;
        }
        
        // 4. 计算Hessian矩阵
        vector<vector<double>> H = hess(x);
        
        // 5. 求解 H * d = -g （使用高斯消元法或直接求逆）
        vector<double> d = solve_linear_system(H, g);
        for (int i = 0; i < n; ++i) {
            d[i] = -d[i];  // d = -H^(-1) * g
        }
        
        // 6. 更新迭代点 X^(k+1) = X^(k) + d^(k)
        for (size_t i = 0; i < x.size(); ++i) {
            x[i] += d[i];
        }
        
        // 打印迭代信息
        cout << iter << "\t";
        for (size_t i = 0; i < x.size(); ++i) {
            printf("%.6lf", x[i]);
            cout << (i == x.size() - 1 ? "" : ", ");
        }
        printf("\t%.6lf\t\t%.6lf\n", func(x), grad_norm);
    }
    
    return x;
}

阻尼牛顿法

引例：标准牛顿法求解

在牛顿法一节里，二次函数一步就能收敛。可实际的目标函数大多不是二次的，这时步长固定为 1 的牛顿法就未必好用了。先看一道例题：

\[f(X) = \frac{1}{4}x_1^4 + \frac{1}{4}x_2^4 - x_1 - x_2,\qquad X^{(0)} = [0.1,\ 0.1]^T,\qquad \varepsilon = 0.01 \]

阻尼牛顿法引例示意图

它的梯度与Hessian矩阵为：

\[\nabla f(X) = \left(\begin{array}{c} x_1^3 - 1 \\ x_2^3 - 1 \end{array}\right),\qquad H(X) = \left(\begin{array}{cc} 3x_1^2 & 0 \\ 0 & 3x_2^2 \end{array}\right) \]

此函数严格凸，极小点唯一：令 \(\nabla f(X) = 0\)，得 \(X^* = [1,\ 1]^T\)，\(f(X^*) = -1.5\)。

而且整条迭代路径上 \(H(X)\) 始终保持正定。

代入 \(X^{(0)}\)：

\[\nabla f(X^{(0)}) = [-0.999,\ -0.999]^T,\qquad H(X^{(0)}) = \left(\begin{array}{cc} 0.03 & 0 \\ 0 & 0.03 \end{array}\right) \]

于是牛顿方向

\[d^{(0)} = -[H(X^{(0)})]^{-1} \nabla f(X^{(0)}) = [33.3,\ 33.3]^T \]

方向没有问题，起始点 \(X^{0}[0.1,\ 0.1]^T\) 与最优解 \([1,\ 1]^T\) 的连线正是 \([1,\ 1]^T\) 方向；问题出在步子上：每个坐标本只需走 \(0.9\) 就能到位，而牛顿法一步迈出去 33.3，足足大了 37 倍：

\[X^{(1)} = X^{(0)} + d^{(0)} = [33.4,\ 33.4]^T,\qquad f(X^{(1)}) = 622170.26 \]

函数值从 \(-0.200\) 直接冲到 \(6.2 \times 10^5\)。运行程序（初始点 0.1 0.1，精度 0.01），完整的迭代过程如下：

步骤 点 X 函数值 f(X) 梯度范数
0 33.400000, 33.400000 622170.256800 1.412799
1 22.266965, 22.266965 122873.087298 52691.764513
2 14.845316, 14.845316 24254.732026 15612.008031
3 9.898390, 9.898390 4780.059295 4625.413538
4 6.602329, 6.602329 936.871811 1370.126348
5 4.409199, 4.409199 180.158593 405.597039
6 2.956612, 2.956612 32.294167 119.811345
7 2.009207, 2.009207 4.129915 35.136683
8 1.422043, 1.422043 -0.799429 10.056459
9 1.112865, 1.112865 -1.458828 2.652586
10 1.011060, 1.011060 -1.499630 0.534923
11 1.000121, 1.000121 -1.500000 0.047443
12 1.000121, 1.000121 -1.500000 0.000511 (收敛)

标准牛顿法先一步“冲上高空”，再花 11 步慢慢滑回谷底，整整 12 次迭代才收敛到 \([1.000121,\ 1.000121]^T\)。

方向是对的，函数是凸的，Hessian是正定的，问题只出在步长——这正是阻尼牛顿法要解决的问题。

问题与改进

上面的例题里，方向没问题、Hessian矩阵也没问题，翻车的唯一原因就是步长。一般地，标准牛顿法存在以下问题：

  1. 当Hessian矩阵不正定时，搜索方向可能不是下降方向
  2. 步长固定为1，可能导致函数值不降反升
  3. 初始点远离极值点时，可能不收敛

阻尼牛顿法（Damped Newton's Method）通过引入一维搜索来确定最优步长，解决了上述问题中的第 2 条，并大幅缓解了第 3 条：只要搜索方向是下降方向，它总能找到合适的步长使函数值下降，迭代过程因此不易失控，对初始点的要求也大大降低。

注意：阻尼牛顿法只调整步长，不改变方向。问题 1（Hessian矩阵不正定导致搜索方向可能不是下降方向）不是调步长能够解决的，需要对方向本身进行修正——详见下一篇文章《3-6 修正牛顿法》。

阻尼牛顿法迭代公式

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

其中：

  • \(d^{(k)} = -[H(X^{(k)})]^{-1} \nabla f(X^{(k)})\) 为牛顿方向
  • \(\alpha^{(k)}\) 通过一维搜索确定，使 \(f(X^{(k)} + \alpha^{(k)} d^{(k)})\) 最小，\(\alpha\) 也称为阻尼因子。

具体步骤

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

  • ① 确定初始点 \(X^{(0)}\)，收敛精度 \(\varepsilon\)，令 \(k=0\)；
  • ② 计算梯度 \(\nabla f(X^{(k)})\) 和 Hessian矩阵 \(H(X^{(k)})\) 以及其逆矩阵 \([H(X^{(k)})]^{-1}\)；
  • ③ 构造搜索方向 \(d^{(k)} = -[H(X^{(k)})]^{-1} \nabla f(X^{(k)})\)；
  • ④ 一维搜索：求最优步长 \(\alpha^{(k)}\)，使 \(f(X^{(k)} + \alpha^{(k)} d^{(k)})\) 最小；
  • ⑤ 更新迭代点：\(X^{(k+1)} = X^{(k)} + \alpha^{(k)} \cdot d^{(k)}\)，\(k=k+1\)；
  • ⑥ 重复②~⑤步，直到：\(\|\nabla f(X^{(k)})\| \le \varepsilon\)。

程序框图

阻尼牛顿法的程序框图

代码示例

阻尼牛顿法 C++ 代码示例： 
// 阻尼牛顿法 求多元函数最小值点
vector<double> damped_newton_method(
    function<double(const vector<double>&)> func,
    function<vector<double>(const vector<double>&)> grad,
    function<vector<vector<double>>>(const vector<double>&)> hess,
    const vector<double>& x0,
    double eps,
    int max_iter = 1000
) {
    vector<double> x = x0;
    int n = x0.size();
    
    cout << "\n迭代过程:" << endl;
    cout << "步骤\t点 X\t\t\t函数值 f(X)\t\t梯度范数\t步长" << endl;
    cout << "--------------------------------------------------------------------" << endl;
    
    for (int iter = 0; iter < max_iter; ++iter) {
        // 1. 计算当前点的梯度
        vector<double> g = grad(x);
        
        // 2. 计算梯度范数
        double grad_norm = norm(g);
        
        // 3. 检查收敛条件
        if (grad_norm < eps) {
            cout << iter << "\t";
            for (size_t i = 0; i < x.size(); ++i) {
                printf("%.6lf", x[i]);
                cout << (i == x.size() - 1 ? "" : ", ");
            }
            printf("\t%.6lf\t\t%.6lf\t\t(收敛)\n", func(x), grad_norm);
            break;
        }
        
        // 4. 计算Hessian矩阵
        vector<vector<double>> H = hess(x);
        
        // 5. 求解牛顿方向 d = -H^(-1) * g
        vector<double> d = solve_linear_system(H, g);
        for (int i = 0; i < n; ++i) {
            d[i] = -d[i];
        }
        
        // 6. 定义一维搜索函数 phi(alpha) = f(x + alpha * d)
        auto phi = [&](double alpha) {
            vector<double> new_x(x.size());
            for (size_t i = 0; i < x.size(); ++i) {
                new_x[i] = x[i] + alpha * d[i];
            }
            return func(new_x);
        };
        
        // 7. 使用进退法+黄金分割法确定最优步长 alpha
        pair<double, double> bracket;
        try {
            bracket = bracket_minimum(phi, 0.0, 0.1);
        } catch (...) {
            printf("[Warning] bracket minimum failed! bracket = [0, 10]\n");
            bracket = make_pair(0.0, 10.0);
        }
        
        double alpha_star = golden_section(phi, bracket.first, bracket.second, eps);
        
        // 8. 更新迭代点 X^(k+1) = X^(k) + alpha* * d^(k)
        for (size_t i = 0; i < x.size(); ++i) {
            x[i] += alpha_star * d[i];
        }
        
        // 打印迭代信息
        cout << iter << "\t";
        for (size_t i = 0; i < x.size(); ++i) {
            printf("%.6lf", x[i]);
            cout << (i == x.size() - 1 ? "" : ", ");
        }
        printf("\t%.6lf\t\t%.6lf\t\t%.6lf\n", func(x), grad_norm, alpha_star);
    }
    
    return x;
}

注：本节（以及前面的牛顿法一节）只列出核心代码片段；包含全部支撑函数的完整可运行程序见下一篇文章《3-6 修正牛顿法》。

回到引例：阻尼牛顿法求解

\(\qquad\)用阻尼牛顿法重解本节开头的例题。牛顿方向仍是 \(d^{(0)} = [33.3,\ 33.3]^T\)，但它不再被全盘接受，而是先沿该方向做一维搜索，求 \(\varphi(\alpha) = f(X^{(0)} + \alpha d^{(0)})\) 的极小点。这条射线恰好穿过最优解 \([1,\ 1]^T\)，对应的理论最优步长是 \(\alpha = 0.9 / 33.3 \approx 0.027\)；程序用进退法 + 黄金分割求得 \(\alpha^{(0)} = 0.025994\)：

\[X^{(1)} = X^{(0)} + 0.025994\, d^{(0)} = [0.965597,\ 0.965597]^T,\qquad f(X^{(1)}) = -1.496530 \]

一步就到了最优解附近，第 2 次迭代再修正一次即可。迭代过程如下：

步骤 点 X 函数值 f(X) 梯度范数 步长
0 0.965597, 0.965597 -1.496530 1.412799 0.025994
1 0.999953, 0.999953 -1.500000 0.140997 0.963872
2 0.999953, 0.999953 -1.500000 0.000201 (收敛)

\(\qquad\)全程函数值单调下降（\(-0.200 \to -1.496530 \to -1.500000\)），没有再出现“过山车”。同样的函数、同样的起点、连方向都一样，标准牛顿法迭代 12 次，阻尼牛顿法 2 次，差别在于步长，这就是阻尼因子的意义。

\(\qquad\)下一篇，我们讲重点探讨一下如果 Hessian 矩阵不正定的情况下，会出现什么问题，以及我们应该如何通过修正来解决出现的问题。

本学习笔记参考资料：

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

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

posted @ 2026-08-01 16:49  H_Elden  阅读(164)  评论(2)    收藏  举报