[最优化技术] 3-1 黄金分割法

3-1 黄金分割法

黄金分割法

简述

黄金分割法是常用的一维搜索试探方法，又称“0.618”法。

黄金分割法的基本思想是：每次缩小后的 新区间长度 与 原区间长度 的比值始终是一个常数，此常数为0.618。也就是说每次的区间缩短率都等于0.618。

• 每次迭代只需计算一个新点的函数值（因为另一个点可复用上一轮的结果），大大减少了函数评估次数。
• 黄金分割比例具有自相似性：每次缩小区间后，新的两点仍保持黄金比例分布。
• 这种对称性和比例不变性避免了数值误差累积，提高了算法的数值稳定性。

证明

由上图所示，为使每次迭代时有一个分割点可以复用，即：上一次的左分割点作为这一次的右分割点；或者上一次的右分割点作为这一次的左分割点，需要满足以下比例成立：

\[1:\lambda = \lambda:(1-\lambda) \]

故：

\[\begin{align*} 1-\lambda &= \lambda^2 \\ \lambda^2 + \lambda - 1 &= 0 \end{align*} \]

取正数解：

\[\lambda = \frac{\sqrt{5}-1}{2} \approx 0.618 \]

此数也就是大名鼎鼎的黄金分割比。

计算公式

\[\begin{align*} a_1 &= b - 0.618(b - a) \\ b_1 &= a + 0.618(b - a) \end{align*} \]

程序框图

提供如下基于黄金分割法进行一维函数搜索的程序框图，注意 图中字母与下文的代码变量并非严格一一对应，仅供参考。

代码示例

对于黄金分割法，给出以下C++函数仅供参考。

// 黄金分割法 求最小值点（横坐标）
double golden_section(function<double(double)> func, double l, double r, double eps) {
	const double ratio = 0.618;
	double a = l, b = r;
	double aa = b - ratio * (b-a);
	double bb = a + ratio * (b-a);
	double y1 = func(aa);
	double y2 = func(bb);
	
	int step = 0;
	printf("[debug] step %d: (%.5lf, %.5lf)\n", step, a, b);
	while (b - a > eps) {			 // 点距准则
		if (y1 >= y2) {
			a = aa;
			aa = bb;                 // 复用上次分割点横坐标
			y1 = y2;                 // 复用上次分割点纵坐标
			bb = a + ratio * (b-a);  // 计算新的分割点
			y2 = func(bb);
		} else {
			b = bb;
			bb = aa;                 // 复用上次分割点横坐标
			y2 = y1;                 // 复用上次分割点纵坐标
			aa = b - ratio * (b-a);  // 计算新的分割点
			y1 = func(aa);
		}
		step++;
		printf("[debug] step %d: (%.5lf, %.5lf)\n", step, a, b);
	}
	return (a + b) / 2;
}

结合上一篇文章所讲的进退法确定初始区间，我们可以轻松求解出任何一维单谷函数的 极小值点 和 极小值。给出如下完整C++代码仅供参考。

C++代码实现：

#include <iostream>
#include <cmath>
#include <functional>
#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)) {
		printf("[debug] left search.\n");
		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");
			printf("[debug] step %d: (%.2lf, %.2lf, %.2lf)\n", iter, x2, x1, x0);
			x0 = x1;
			x1 = x2;
			h *= 2;
			x2 = x1 - h; // 步长加倍
		}
		printf("[debug] step %d: (%.2lf, %.2lf, %.2lf)\n", iter + 1, x2, x1, x0);
		return make_pair(x2, x0);
	}
	// 递减，向右搜索
	else {
		printf("[debug] right search.\n");
		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");
			printf("[debug] step %d: (%.2lf, %.2lf, %.2lf)\n", iter, x0, x1, x2);
			x0 = x1;
			x1 = x2;
			h *= 2;
			x2 = x1 + h; // 步长加倍
		}
		printf("[debug] step %d: (%.2lf, %.2lf, %.2lf)\n", iter + 1, x0, x1, x2);
		return make_pair(x0, x2);
	}
}


// 黄金分割法 求最小值点（横坐标）
double golden_section(function<double(double)> func, double l, double r, double eps) {
	const double ratio = 0.618;
	double a = l, b = r;
	double aa = b - ratio * (b-a);
	double bb = a + ratio * (b-a);
	double y1 = func(aa);
	double y2 = func(bb);
	
	int step = 0;
	printf("[debug] step %d: (%.5lf, %.5lf)\n", step, a, b);
	while (b - a > eps) {			// 点距准则
		if (y1 >= y2) {
			a = aa;
			aa = bb;
			y1 = y2;
			bb = a + ratio * (b-a);
			y2 = func(bb);
		} else {
			b = bb;
			bb = aa;
			y2 = y1;
			aa = b - ratio * (b-a);
			y1 = func(aa);
		}
		step++;
		printf("[debug] step %d: (%.5lf, %.5lf)\n", step, a, b);
	}
	return (a + b) / 2;
}

int main() {
	// 定义一个一元凸函数
	auto f = [](double x) {
		return x * x - 2 * x - 5;
	};

	double a0, h, eps;
	cout << "输入初始点a0: ";
	cin >> a0;
	cout << "输入初始步长h: ";
	cin >> h;
	cout << "输入误差精度epsilon: ";
	cin >> eps;

	pair<double, double> range;

	try {
		range = bracket_minimum(f, a0, h);
		cout << "初始搜索区间为：(" << range.first << ", " << range.second << ").\n";
	} catch (const exception& e) {
		cerr << "Error: " << e.what() << endl;
		return 1;
	}

	double res = golden_section(f, range.first, range.second, eps);
	cout << "最小值点 为：" << res << " ± " << eps / 2 << endl;
	cout << "最小值   为：" << f(res);
	return 0;
}

扩展

此外，类似黄金分割法的等比例缩小区间的方法还有很多。

代码结构如下所示。

C++代码示例：

double ratio_section(function<double(double)> func, double l, double r, double eps) {
    const double ratio = 1.0/3;
	double a = l, b = r, aa, bb;
	int step = 0;
	printf("[debug] step %d: (%.5lf, %.5lf)\n", step, a, b);
	while (b - a > eps) {			// 点距准则
		aa = a + ratio * (b - a);	// 确定左分割点
		bb = b - ratio * (b - a);	// 确定右分割点
		step++;
		if (func(aa) < func(bb))
			b = bb;
		else
			a = aa;
		printf("[debug] step %d: (%.5lf, %.5lf)\n", step, a, b);
	}
	return (a + b) / 2;
}

当然，函数中的 ratio 可以改成任何你喜欢的数字，当它为 \(\frac{1}{3}\) 时，例如上面代码所示，这就变成了算法竞赛中经典的 三分法。

如果把 ratio 设置为 1 - 0.618 那么很显然，这就又变成了 黄金分割法 了。当然，它并没有利用好黄金分割法 每次迭代只新增一个分割点 的优势，这个代码是需要计算两次的。

---

本学习笔记参考资料：

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

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