[最优化技术] 3-5 牛顿法与阻尼牛顿法
3-5 牛顿法与阻尼牛顿法
牛顿法
简述
\(\qquad\)牛顿法(Newton's Method)是一种利用目标函数的二阶导数信息(Hessian矩阵)进行优化的方法。与最速下降法仅使用一阶梯度信息不同,牛顿法通过二次泰勒展开近似原函数,能够更准确地预测极值点的位置,因此具有更快的收敛速度(二阶收敛)。
\(\qquad\)牛顿法的基本思想是:在当前点处用二次函数近似原目标函数,然后直接求出该二次函数的极小值点作为下一次迭代点。
原理
对于多元函数 \(f(X)\),在点 \(X^{(k)}\) 处进行二阶泰勒展开:
其中 \(H(X^{(k)})\) 为 Hessian矩阵(海森矩阵):
注意:当函数的所有二阶混合偏导数连续时,有 \(\frac{\partial^2 f}{\partial x_i \partial x_j} = \frac{\partial^2 f}{\partial x_j \partial x_i}\),即 H矩阵对称。对于多项式函数,Hessian矩阵必然对称。
对于二次函数求极值,令泰勒展开的一阶导数为零:
解得牛顿法迭代公式:
Hessian矩阵的性质
Hessian矩阵的正定性决定了极值点的性质:
- 正定(所有特征值 > 0)→ 局部极小点
- 负定(所有特征值 < 0)→ 局部极大点
- 不定(有正有负特征值)→ 鞍点
- 半定(有零特征值)→ 需要高阶检验
对于二元二次函数,Hessian矩阵为:
二阶矩阵求逆公式:对于矩阵 \(A = \left(\begin{array}{cc} a & b \\ c & d \end{array}\right)\),当 \(ad - bc \neq 0\) 时:
牛顿法分析
特点:
- 原函数为二次函数:一步求出极值点。
- 原函数为非二次函数:先通过二阶泰勒展开变为二次函数,因此求出的为极值点的近似值。
优点:
- 收敛速度最快(二阶收敛)
- 对于二次函数,一次迭代即可收敛
缺点:
- Hessian矩阵计算困难(需要计算所有二阶偏导数)
- Hessian矩阵求逆计算量大
- 当Hessian矩阵不正定时,可能不收敛或收敛到鞍点
具体步骤
对于给定初始点 \(X^{(0)}\) 和收敛精度 \(\varepsilon\):
- ① 确定初始点 \(X^{(0)}\),收敛精度 \(\varepsilon\),令 \(k=0\);
- ② 计算梯度 \(\nabla f(X^{(k)})\) 和 Hessian矩阵 \(H(X^{(k)})\);
- ③ 判断Hessian矩阵是否正定。若不正定,可采用修正策略(见阻尼牛顿法);
- ④ 构造搜索方向 \(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)\) 的梯度函数和Hessian矩阵:
梯度:
Hessian矩阵:
注意:对于二次函数,Hessian矩阵为常数矩阵,不随 \(X\) 变化。
代入 \(X^{(0)} = [1,1]^T\),得:
计算梯度的模长:
不满足收敛条件,迭代继续。
第 1 次迭代 (\(k=0\)):
- 计算Hessian矩阵的逆:
行列式:\(\det(H) = 2 \times 4 - (-2) \times (-2) = 8 - 4 = 4\)
- 计算牛顿方向:
- 更新迭代点:
- 检验收敛条件:
计算 \(X^{(1)}\) 处的梯度:
梯度的模长:
满足收敛条件,迭代结束!
最终得到最优解:\(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;
}
阻尼牛顿法
问题与改进
标准牛顿法存在以下问题:
- 当Hessian矩阵不正定时,搜索方向可能不是下降方向
- 步长固定为1,可能导致函数值不降反升
- 初始点远离极值点时,可能不收敛
阻尼牛顿法(Damped Newton's Method)通过引入一维搜索来确定最优步长,解决了上述问题。
阻尼牛顿法迭代公式
其中:
- \(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;
}
完整示例
C++代码实现(完整):
#include <iostream>
#include <vector>
#include <cmath>
#include <functional>
#include <iomanip>
#include <utility>
using namespace std;
// 计算向量范数
double norm(const vector<double>& v) {
double sum = 0.0;
for (double vi : v) {
sum += vi * vi;
}
return sqrt(sum);
}
// 进退法求初始区间
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 in left expansion");
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 in right expansion");
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;
}
// 求解线性方程组 H * x = b (高斯消元法)
vector<double> solve_linear_system(vector<vector<double>> H, vector<double> b) {
int n = b.size();
// 增广矩阵
for (int i = 0; i < n; ++i) {
H[i].push_back(b[i]);
}
// 高斯消元
for (int i = 0; i < n; ++i) {
// 找主元
int max_row = i;
for (int k = i + 1; k < n; ++k) {
if (abs(H[k][i]) > abs(H[max_row][i])) {
max_row = k;
}
}
// 交换行
swap(H[i], H[max_row]);
// 消元
for (int k = i + 1; k < n; ++k) {
double factor = H[k][i] / H[i][i];
for (int j = i; j <= n; ++j) {
H[k][j] -= factor * H[i][j];
}
}
}
// 回代
vector<double> x(n);
for (int i = n - 1; i >= 0; --i) {
x[i] = H[i][n];
for (int j = i + 1; j < n; ++j) {
x[i] -= H[i][j] * x[j];
}
x[i] /= H[i][i];
}
return x;
}
// 阻尼牛顿法
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) {
vector<double> g = grad(x);
double grad_norm = norm(g);
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;
}
vector<vector<double>> H = hess(x);
vector<double> d = solve_linear_system(H, g);
for (int i = 0; i < n; ++i) {
d[i] = -d[i];
}
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);
};
pair<double, double> bracket;
try {
bracket = bracket_minimum(phi, 0.0, 0.1);
} catch (...) {
printf("[Warning] bracket minimum failed!\n");
bracket = make_pair(0.0, 10.0);
}
double alpha_star = golden_section(phi, bracket.first, bracket.second, eps);
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;
}
// ==========================================================
// 示例:二元二次凸函数
// f(X) = x1^2 + 2*x2^2 - 4*x1 - 2*x1*x2
// ==========================================================
int main() {
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 = [](vector<double> x) {
return x[0] * x[0] + 2 * x[1] * x[1] - 4 * x[0] - 2 * x[0] * x[1];
};
// 梯度函数
auto grad = [](vector<double> x) {
vector<double> g(2);
g[0] = 2 * x[0] - 2 * x[1] - 4;
g[1] = 4 * x[1] - 2 * x[0];
return g;
};
// Hessian矩阵函数
auto hess = [](vector<double> x) {
vector<vector<double>> H(2, vector<double>(2));
H[0][0] = 2; // ∂²f/∂x1²
H[0][1] = -2; // ∂²f/∂x1∂x2
H[1][0] = -2; // ∂²f/∂x2∂x1
H[1][1] = 4; // ∂²f/∂x2²
return H;
};
double x1_init, x2_init, eps;
cout << "输入初始点 x1: ";
cin >> x1_init;
cout << "输入初始点 x2: ";
cin >> x2_init;
cout << "输入精度 eps: ";
cin >> eps;
vector<double> x0 = {x1_init, x2_init};
vector<double> result = damped_newton_method(f, grad, hess, x0, eps);
cout << "\n最终结果:" << endl;
cout << "最优解 X* = [";
for (size_t i = 0; i < result.size(); ++i) {
cout << result[i] << (i == result.size() - 1 ? "" : ", ");
}
cout << "]^T" << endl;
cout << "目标函数值 f(X*) = " << f(result) << endl;
vector<double> final_grad = grad(result);
double final_norm = norm(final_grad);
cout << "最终梯度范数 = " << final_norm << endl;
return 0;
}
总结
牛顿法与阻尼牛顿法的比较:
| 特性 | 牛顿法 | 阻尼牛顿法 |
|---|---|---|
| 步长 | 固定为1 | 通过一维搜索确定 |
| 收敛性 | 局部收敛,要求初始点靠近极值点 | 全局收敛性更好 |
| 计算量 | 较小(无需一维搜索) | 较大(需一维搜索) |
| 适用场景 | 二次函数或初始点较好时 | 一般情况,更稳健 |
牛顿法的优缺点:
优点:
- 收敛速度最快(二阶收敛)
- 对于二次函数,一次迭代即可收敛到精确解
- 充分利用了函数的二阶信息
缺点:
- 需要计算Hessian矩阵(所有二阶偏导数)
- 需要求解线性方程组或矩阵求逆,计算量大
- 当Hessian矩阵不正定或奇异时,可能不收敛
- 对初始点要求较高(标准牛顿法)
阻尼牛顿法的改进:
- 通过一维搜索确定最优步长,保证了函数值下降
- 具有更好的全局收敛性
- 即使初始点远离极值点,也能稳定收敛
适用场景:
- 目标函数的二阶导数容易计算
- 对收敛速度要求较高
- 问题规模适中(Hessian矩阵不太大)
\(\qquad\)一般地,将牛顿法和阻尼牛顿法统称为牛顿型方法。牛顿型方法总体上迭代次数较少、计算速度较快。但是这类方法的主要缺点是每次迭代都要计算函数的二阶导数矩阵,并对该矩阵求逆。当目标函数的维数高时,算量和存储量大的缺点尤为明显。最速下降法的收敛速度比牛顿型方法慢,而牛顿型方法又存在上述缺点。
\(\qquad\)至此,我们已经介绍了两种基于梯度信息的优化方法:最速下降法原理简单但收敛较慢,牛顿法收敛速度快但需要计算Hessian矩阵。在接下来的文章中,我们将首先介绍坐标轮换法——一种无需计算梯度的直接搜索方法,它通过沿坐标轴方向依次搜索来寻找极值点,实现更为简单。之后,我们将介绍共轭梯度法,它巧妙地结合了最速下降法和牛顿法的优点,既避免了Hessian矩阵的计算,又克服了最速下降法的锯齿现象,具有较快的收敛速度。
本学习笔记参考资料:
[1] 白清顺, 孙靖民, 梁迎春. 机械优化设计 第7版[M]. 北京: 机械工业出版社, 2024. ISBN: 978-7-111-75103-8
[2] 武汉理工大学《最优化技术B》课程课件,授课教师:颜彬老师

本文详细介绍多元函数无约束优化中的牛顿法与阻尼牛顿法。该方法利用目标函数的二阶导数信息(Hessian矩阵)构造二次近似,通过求解线性方程组确定搜索方向,具有二阶收敛速度。文章系统阐述了Hessian矩阵的性质、牛顿法的迭代公式及其几何意义,并通过手算实例展示了二次函数一步收敛的特性。针对标准牛顿法的局限性,引入阻尼牛顿法通过一维搜索确定最优步长,提高了算法的全局收敛性与稳定性。
浙公网安备 33010602011771号