分类搜索方法与导数方法——《数值计算方法》

Dropdrag

已于 2022-01-25 22:00:02 修改

阅读量1k

点赞数 1

文章标签：分类机器学习算法

于 2022-01-25 20:04:23 首次发布

本文链接：https://blog.csdn.net/qq_52867930/article/details/122690907

版权

《数值计算方法》系列总目录

第一章误差序列实验
第二章非线性方程f(x)=0求根的数值方法
第三章 CAD模型旋转和AX=B的数值方法
第四章插值与多项式逼近的数值计算方法
第五章曲线拟合的数值方法
第六章数值微分计算方法
第七章数值积分计算方法
第八章数值优化方法

一、算法原理

对于求解一元函数的极小值点，常用的方法即为导数法。即根据极值点的性质，在极值点处函数的导数为 $0$ 。因此，我们给出极值点的定义。如果存在包含 $p$ 的开区间 $I$ ，使得对所有 $x \in I$ ，有 $f (p) \leq f (x)$ ，则称函数 $f$ 在 $x = p$ 处有局部极小值。类似地，如果对所有 $x \in I$ ，有 $f (x) \leq f (p)$ ，则称 $f$ 在 $x = p$ 处有局部极大值。如果 $f$ 在点 $x = p$ 处有局部极大值或极小值，则称 $f$ 在 $x = p$ 处有局部极值。因此，我们可以得出，对于极小值来说，在极小值点的一个领域内，极小值点左侧的点的导数小于 $0$ ，在右侧则大于 $0$ ，这是由函数的导数以及单调性可知。对于极大值点来说则正好相反。因此我们可以得到通过一阶导数来判断函数极小值点的方法。既有如下定理：设 $f (x)$ 在区间 $I = [a, b]$ 上连续，并设除 $x = p$ 处外，对所有 $x \in (a, b)$ 都有定义。则有
i. 若在 $(a, p)$ 上 $f^{'} (x) < 0$ ，而在 $(p, b)$ 上 $f^{'} (x) > 0$ ，则 $f (p)$ 是局部极小值。
ii. 若在 $(a, p)$ 上 $f^{'} (x) > 0$ ，而在 $(p, b)$ 上 $f^{'} (x) < 0$ ，则 $f (p)$ 是局部极大值。
除此之外，根据函数二阶导数的定义以及一阶导数的单调性可知，还可以通过二阶导数来判断函数的极值点。因此有如下定理：设f在区间 $I = [a, b]$ 上连续，并且 $f^{'} (x)$ 和 $f^{''} (x)$ 区间 $(a, b)$ 上有定义，又设 $p \in (a, b)$ 是关键点，即 $f^{'} (x) = 0$ 。则有
i. 若 $f^{''} (x) > 0$ ，则 $f (p)$ 是f的一个局部极小值。
ii. 若 $f^{''} (x) < 0$ ，则 $f (p)$ 是f的一个局部极大值。
iii. 若 $f^{''} (x) = 0$ ，则结果不确定。
因此，可以通过函数的导数来求解函数的极值点。然而是否可以通过在函数导数未知的情况下求解函数的极值点？分类搜索法可以在函数导数未知的情况下求解函数的极值点，即通过极值点的定义，在一个已知含有函数极值点的区间上通过迭代以及极值点的判定依据使区间收敛，从而得到在误差精度内的极值点。这种方法叫做分类搜索法。
然而如何选取区间上的迭代点以及搜索方式能够使得算法效率更高是一个值得考虑的问题。要减少函数求值的次数，确定在哪里求 $f (p)$ 值的好策略非常重要。其中黄金分割搜索法和斐波那契搜索法是两种有效的方法。使用这两种方法来求 $f (p)$ 的极小值必须满足特定的条件，以保证在给定的区间内有合适的极小值。因此在搜索之前我们还应该直到一个区间内是否存在极值点，则可以通过函数单峰的定义得到，即如果存在惟一的 $p \in I$ ，使得（1） $f (p)$ 在 $[a, p]$ 上递减(2) $f (p)$ 在 $[p, b]$ 上递增，则函数 $f (p)$ 在 $I = [a, b]$ 上是单峰的。因此我们可以在单峰区间对函数进行分类搜索法。

1、黄金分割搜索法

如果已知 $f (x)$ 在 $[a, b]$ 上是单峰的，则有可能找到该区间的一个子区间， $f (x)$ 在该子区间上取得极小值。一种方法是选择两个内点 $c < d$ ，这样就有 $a < c < d < b$ 。f（x）的单峰特性保证了函数值 $f (c)$ 和 $f (d)$ 小于 $max\{f(a)，f(b)\}$ 。这时有两种情况要考虑，即

如果 $f (c) \leq f (d)$ ，则最小值必然落在子区间 $[a, d]$ 中。可将 $b$ 替换为 $d$ ，并在新区间 $[a, d]$ 上继续搜索。如果 $f (d) < f (c)$ ，则最小值必然落在子区间 $[c, b]$ 中。将 $a$ 替换为 $c$ ，并在新区间 $[c, b]$ 中继续搜索。选择内点 $c$ 和 $d$ ，使得区间 $[a, c]$ 与 $[d, b]$ 对称，即 $b - d = c - a$ ，其中
$r)b\\ d = b - (1 - r)(b - a) = (1 - r)a + rb$ 并且 $1 / 2 < r < 1$ （以保证 $c < d$ ）。希望 $r$ 在每个子区间上保持为常数。并且，旧的内点中的一个应该是新子区间的一个内点。而另一个则应该是新子区间的一个端点。这样，在每次迭代中只需要找一个新的点，只需要一次新的函数求值计算。这样一来，搜索迭代的速度就会大幅度提高。如果 $f (c) \leq f (d)$ ，并且只能进行一次新的函数求值，则必须要求
$\frac{{d - a}}{{b - a}} = \frac{{c - a}}{{d - a}}$ 通过计算解得 $\frac{{ - 1 \pm \sqrt 5 }}{2}$ （黄金分割比例）。类似地，如果 $f (d) < f (c)$ ，则解得 $\frac{{ - 1 \pm \sqrt 5 }}{2}$ 。这样一来，每次搜索迭代过程在判定搜索方向后，就可以只计算一次新的迭代点即可，这样一来就大大提高了计算的效率。其中，每次搜索所得区间的比例应为黄金分割比。

2、斐波那契搜索法

在黄金分割搜索法中，第一次迭代中进行了两次函数求值，而在后续的每次迭代中则只进行一次函数求值。 $r$ 的值对每个子区间相同，当 ${b_k} - {a_k}|$ 或 $f({b_k}) - f({a_k})|$ 满足预定义的容差时，迭代结束。斐波那契搜索法与黄金分割搜索法的不同之处在于， $r$ 的值并不是常数。并且，子区间数（迭代数）是由指定的容差决定的。斐波那契搜索基于公式
${F_0} = 0,{F_1} = 1,{F_n} = {F_{n - 1}} + {F_{n - 2}}$ $n = 2 ， 3 ， \dots$ 定义的斐波那契序列 $\{ {F_k}\} _0^\infty$ 。因此斐波那契数为 $0 ， 1 ， 1 ， 2 ， 3 ， 5 ， 8 ， 13 ， 21 ， \dots$ 。设给定函数 $f (x)$ 在区间 ${a_0},{b_0}]$ 上是单峰的。与黄金分割搜索法中一样，也要选择一个值 $r$ ，其中 $1 / 2 < r < 1$ ，与黄金分割比搜索法思想一样，要使得内点 ${c_0},{d_0}$ 可以在下一个子区间上使用，从而只需一次新的函数求值计算。不失一般性，设 $f({c_0}) > f({d_0})$ ，它满足 ${a_1} = {a_0},{b_1} = {d_0},{d_1} = {c_0}$ 。如果只能进行一次新的函数求值计算，则为子区间 $[a, b]$ 选择 $r$ 使得
${d_0} - {c_0} = {b_1} - {d_1}\\ (2{r_0} - 1)({b_0} - {a_0}) = (1 - {r_1})({b_1} - {a_1})\\ (2{r_0} - 1)({b_0} - {a_0}) = (1 - {r_1}){r_0}({b_0} - {a_0})\\ (2{r_0} - 1) = (1 - {r_1}){r_0}\\ \Rightarrow {r_1} = \frac{{1 - {r_0}}}{{{r_0}}}$ 代入公式，可以得到
${r_1} = \frac{{1 - \frac{{{F_{n - 1}}}}{{{F_n}}}}}{{\frac{{{F_{n - 1}}}}{{{F_n}}}}} = \frac{{{F_n} - {F_{n - 1}}}}{{{F_{n - 1}}}} = \frac{{{F_{n - 2}}}}{{{F_{n - 1}}}}$ 因此有，斐波那契搜索可以从开始，用
${r_k} = \frac{{{F_{n - 1 - k}}}}{{{F_{n - k}}}},k = 1,2,...,n - 3$ 注意 $r_3=F_2/F_3=1/2$ ，因此这一步中无需增加新的点。因此，这一过程总共需要 $(n - 3) ＋ 1 = n - 2$ 步。将第 $k$ 个子区间的长度按因子 ${r_k} = \frac{{{F_{n - 1 - k}}}}{{{F_{n - k}}}}$ 缩减，得到第 $(k + 1)$ 个子区间。最后一个子区间的长度为
$\frac{{{F_{n - 1}}{F_{n - 2}}...{F_2}}}{{{F_n}{F_{n - 1}}...{F_3}}}({b_0} - {a_0}) = \frac{{{F_2}}}{{{F_n}}}({b_0} - {a_0}) = \frac{1}{{{F_n}}}({b_0} - {a_0})$ 对于最后一个区间，通过设置容差即可求出 ${F_n}$ 的大小，通过斐波那契数列的递推公式即可确定 $n$ 即迭代次数，因此斐波那契搜索法的迭代次数由设置的初始容差决定。

3、利用导数求解极小值

利用导数求解极小值点的核心思想就是利用基于函数在一个包含极小值点的区间的三个关键点构造的拉格朗日多项式逼近函数，并求解出多项式的导数以及导数取得0时的近似极小值点表达式，通过迭代使得近似极小值点表达式逐渐向真实极小值点收敛，从而得到极小值点。因此在每次迭代过程中就需要确定基于哪三个点构造拉格朗日多项式。
因此在一个区间内，首先搜索满足条件的区间的初始点以及步长。设有三个点， ${p_0},{\rm{ }}{p_1} = {p_0} + h,{\rm{ }}{p_2} = {p_0} + 2h$ ，满足
$f({p_0}) > f({p_1}),{\rm{ }}f({p_1}) < f({p_2})$ 则区间内包含极小值点。设 $f'({P_0})$ <0，则 $P_0<P$ ，且应该选择正的步长 $h$ 。容易找到 $h$ ，使3个点满足公式。在三个点中，如果满足条件 $a + 1 < b$ ，则令 $h = 1$ ;如果条件不满足，则令 $h = 1 / 2$ ，依次类推。
情况（i）∶若式满足则结束。
情况（ii）∶若 $f({P_0})$ > $f({P_1})$ 且 $f({P_1})$ > $f({P_2})$ ，则说明 $P_2$ < $P$ 。那么，需要检测更靠右的点。将步长加倍，并重复检测过程。
情况（iii）∶若 $f({P_0})≤f({P_1})$ ，表明 $h$ 太大，已经跳过了 $p$ 。那么，需要检测更靠近的点。将步长减半，并重复检测过程。
其中当 $f({P_0})>0$ 时，应该选择负的步长 $h$ ，并用类似于上面的分类进行处理。通过以上步骤，我们即得到满足条件的三个点，因此我们可以构造拉格朗日多项式。既有如下函数近似表达式
$\frac{{{y_0}(x - {p_1})(x - {p_2})}}{{2{h^2}}} - \frac{{{y_1}(x - {p_0})(x - {p_2})}}{{{h^2}}} + \frac{{{y_2}(x - {p_0})(x - {p_1})}}{{2{h^2}}}$ 通过对近似表达式求解一阶导数，并将满足条件的和步长 $h$ 带入表达式，可以得到
$\frac{{{y_0}(2x - {p_1} - {p_2})}}{{2{h^2}}} - \frac{{{y_1}(2x - {p_0} - {p_2})}}{{{h^2}}} + \frac{{{y_2}(2x - {p_0} - {p_1})}}{{2{h^2}}}\\ 0 = \frac{{{y_0}(2({p_0} + {h_{\min }}) - {p_1} - {p_2})}}{{2{h^2}}} - \frac{{{y_1}(2({p_0} + {h_{\min }}) - {p_0} - {p_2})}}{{{h^2}}} + \frac{{{y_2}(2({p_0} + {h_{\min }}) - {p_0} - {p_1})}}{{2{h^2}}}$ 最后得到极小值点表达式，
${h_{\min }} = \frac{{h(4{y_1} - 3{y_0} - {y_2})}}{{4{y_1} - 2{y_0} - 2{y_2}}}$ 因此可以得到步长表达式，因此可以知道，是比点更接近极小值点的点。因此通过重复以上步骤迭代，可以得到一系列收敛于极小值点的迭代点序列，最终在误差精度范围内停止迭代，可以求解出最后的极小值点精确值。

二、实验内容及核心算法代码

1、分类搜索方法与导数方法原理实现

由第一节原理部分可知，分类搜索方法和导数方法最重要的区别在于前者适用范围更广，不要求函数的可导性，然而分类搜索方法往往搜索速度更慢，即需要更多的迭代次数才能求解极小值点。两种搜索方法的原理核心一致，不同的是其每次选取的区间缩小比例不同，因此其泛式流程如下：

选取包含极小值点的区间 ${a_k},{b_k}]$ ，求解区间端点的函数值。
确定第 $k$ 次迭代的区间缩小比例，即求解第 $k$ 次迭代点的比例系数 $r$ （依照两种方法不同的选取原则选取）。
根据系数r求解新的迭代点 ${c_k},{d_k}({c_k} < {d_k})$ 。其中若满足 $f({c_k}) < f({d_k})$ ，则将区间修改为 ${a_k},{d_k}]$ ，令下一个迭代点的 ${d_{k + 1}} = {c_k}$ ；若满足 $f({c_k}) > f({d_k})$ ，则将区间修改为 $f({c_k}) > f({d_k})$ ，令下一个迭代点的 ${c_{k + 1}} = {d_k}$ 。
求解在每个迭代点的函数值，比较是满足精度条件。当满足时停止迭代；否则重复步骤1-3。

对于导数法求解极小值点的方法则更为复杂一点，因为每次迭代过程还涉及极限值点区间以及步长的搜索。其流程大致为：

选取搜索区间的初始值点 $P_0$ 和步长 $h$ ，计算 ${P_0},{P_0} + h,{P_0} + 2h$ 点的函数值。
若 $f({p_0}) > f({p_1}),{\rm{ }}f({p_1}) < f({p_2})$ 满足则结束；若 $f({P_0})>f({P_1})$ 且 $f({P_1})>f({P_2})$ ，将 $h = 2 h$ ，并重复检测过程；若 $f({P_0})<f({P_1})$ ， $h = h / 2$ ，并重复检测过程.
找到初始值点 $P_0$ 和步长h后，利用公式计算第 $k$ 次搜索的优化步长 ${h_{\min }}$ 。
判断 ${h_{\min }}$ 是否满足精度要求，若满足则停止；若不满足，则继续重复步骤1-3，并令 ${P_0} = {P_0} + {h_{\min }}$ 。

通过以上流程图，我们选取以下示例对三个方法进行测试。

利用黄金分割比搜索法、斐波那契搜索法和二次导数逼近法求解以下函数的极小值点：
$(a){\rm{ }}f(x) = {e^x} + 2x + \frac{{{x^2}}}{2};[ - 2.4, - 1.6]\\ (b){\rm{ }}f(x) = - \sin (x) - x + \frac{{{x^2}}}{2};[0.8,1.6]\\ (c){\rm{ }}f(x) = \frac{{{x^2}}}{2} - 4x - x\cos (x);[0.5,2.5]\\ (d){\rm{ }}f(x) = {x^3} - 5{x^2} + 23;[1,5]$ 由以上算法步骤可以知道，对于给定区间上黄金分割比搜索法不需要做任何处理即可直接进行迭代；而斐波那契搜索法则需要先根据确定的误差计算迭代的次数；二次逼近搜索法则需要计算先判断区间端点的一阶导数的正负。
其核心函数算法代码为：

//---------------------------------Golden_section--------------------------------
float GOLDEN_RATIO = (-1. + sqrt(5)) / 2.;
void CmpExtremumbyGolden_section(float (*fun)(float), float a0, float b0, float** funval, const float absci_tol, const float ordin_tol, int& t_clock, const int Max_Iter)
{
	//set the initial ck,dk
	float ck = GOLDEN_RATIO * (a0 - b0) + b0, dk = a0 + GOLDEN_RATIO * (b0 - a0);
	//the f(ck),f(dk)
	float y_ck, y_dk;
	do
	{
		if (t_clock == Max_Iter)
		{
			cout << "Can not find the extremum point within the Max_Iter!" << endl;
			return;
		}
		y_ck = fun(ck);
		y_dk = fun(dk);
		funval[t_clock][0] = t_clock;
		funval[t_clock][1] = a0;
		funval[t_clock][2] = ck;
		funval[t_clock][3] = dk;
		funval[t_clock][4] = b0;
		funval[t_clock][5] = y_ck;
		funval[t_clock][6] = y_dk;
		//judge which side to search
		if (y_ck < y_dk)
		{
			b0 = dk;
			dk = ck;
			ck = GOLDEN_RATIO * (a0 - b0) + b0;
		}
		else
		{
			a0 = ck;
			ck = dk;
			dk = a0 + GOLDEN_RATIO * (b0 - a0);
		}
		++t_clock;
	} while (abs(a0 - b0) > absci_tol || abs(y_ck - y_dk) > ordin_tol);
}
//---------------------------------Fibonacci--------------------------------
//cmp the Fibonacci sequence length by the eps,and return the sequence
float* CmpFibonaciiLength(const float a0, const float b0, int& length, const float eps)
{
	length = 2;
	float Fn = (b0 - a0) / eps;
	float fk_1 = 0, fk = 1, t_fk;
	float* FiboSeque;
	do
	{
		t_fk = fk + fk_1;
		fk_1 = fk;
		fk = t_fk;
		++length;
	} while (t_fk <= Fn);
	FiboSeque = MallocArr(length);
	FiboSeque[0] = 0;
	FiboSeque[1] = 1;
	int t_clock = 2;
	fk_1 = 0, fk = 1;
	do
	{
		t_fk = fk + fk_1;
		fk_1 = fk;
		fk = t_fk;
		FiboSeque[t_clock] = t_fk;
		++t_clock;
	} while (t_clock < length);
	return FiboSeque;
}
void CmpExtremumbyFibonacci(float (*fun)(float), float a0, float b0, float** funval, float* fibo_seque, const int length, const float e)
{
	int iter = length - 1;
	int t_clock = 0;
	float y_ck, y_dk;
	float fibo_ratio = fibo_seque[iter - t_clock - 1] / fibo_seque[iter - t_clock];
	//initial the fibo-ratio
	float ck = a0 * fibo_ratio + (1 - fibo_ratio) * b0, dk = fibo_ratio * b0 + (1 - fibo_ratio) * a0;
	bool done = false;
	do
	{
		//note the value
		funval[t_clock][0] = t_clock;
		funval[t_clock][1] = a0;
		funval[t_clock][2] = ck;
		funval[t_clock][3] = dk;
		funval[t_clock][4] = b0;
		++t_clock;
		float fibo_ratio = fibo_seque[iter - t_clock - 1] / fibo_seque[iter - t_clock];
		if (iter - t_clock == 3)
		{
			done = true;
		}
		y_ck = fun(ck);
		y_dk = fun(dk);
		//determine which side to divide
		if (y_ck < y_dk)
		{
			b0 = dk;
			dk = ck;
			if (done)
			{
				fibo_ratio -= e;
			}
			ck = ck = a0 * fibo_ratio + (1 - fibo_ratio) * b0;
		}
		else
		{
			a0 = ck;
			ck = dk;
			if (done)
			{
				fibo_ratio += e;
			}
			dk = fibo_ratio * b0 + (1 - fibo_ratio) * a0;
		}
	} while (iter - t_clock >= 3);
}
//---------------------------------SecAppro--------------------------------
//cmp the initial step-length based on the interval length
float CmpinitialH(float a0, float b0)
{
	if (a0 + 1 < b0)
	{
		return 1.;
	}
	else
	{
		return .5;
	}
}
//cmp the p0 and h which meet the condition interval including extremum
float CmptestvalatExtre(float (*fun)(float), float& p0, const float h0)
{
	float h = h0;
	float p1, p2, dy_p0, dy_p1, dy_p2;
	do
	{
		p1 = p0 + h;
		p2 = p1 + h;
		dy_p0 = fun(p0);
		dy_p1 = fun(p1);
		dy_p2 = fun(p2);
		//determine the step-length of next iteration
		if (dy_p0 >= dy_p1 && dy_p1 <= dy_p2)
		{
			break;
		}
		else if (dy_p0 > dy_p1 && dy_p1 > dy_p2)
		{
			h *= 2.;
		}
		else
		{
			h /= 2.;
		}
	} while (true);
	return h;
}
//--------------Golden_ratio_search------------------
void Golden_ratio_search(float (*fun)(float), string filename, float a0, float b0)
{
	int t_clock = 0;
	const int Max_Iter = 50;
	const float absci_tol = 1e-6;
	const float ordin_tol = 1e-8;
	float** funval;
	funval = MallocArr(Max_Iter, 7);
	//cmp the Exatremum point
	CmpExtremumbyGolden_section(fun, a0, b0, funval, absci_tol, ordin_tol, t_clock, Max_Iter);
	//save the data
	string title[] = { "k","ak","ck","dk","bk","f(ck)","f(dk)" };
	SaveTabletoFile(filename, funval, title, t_clock, 7);
}
//--------------Fibonacci_search------------------
void Fibonacci_search(float (*fun)(float), string filename, float a0, float b0)
{
	int length = 0;
	const float eps = 1e-4, e = 1e-2;
	float* FiboSeque;
	float** funval;
	FiboSeque = CmpFibonaciiLength(a0, b0, length, eps);
	funval = MallocArr(length, 5);
	CmpExtremumbyFibonacci(fun, a0, b0, funval, FiboSeque, length, e);
	string title[] = { "k","ak","ck","dk","bk" };
	SaveTabletoFile(filename, funval, title, length - 3, 5);
}

四、结论

由以上实验数据可以发现，通常情况下黄金分割比方法是搜索最慢的一种方法。相较于同类型方法的斐波那契搜索法，在相同误差精度要求下，黄金分割比方法的迭代次数要更多，但是二者结果一致。然而，其中收敛最快的是二次逼近方法，其原因是二次逼近方法是基于拉格朗日多项式推导出来的，其搜索方向非常准确，并且在每次搜索前还有确定搜索区间的步骤，因此去搜索范围精度大大提升，使得收敛速度更快。

声明：本系列《数值计算方法》所有章节源代码均为作者编写，大家可以自行下载，仅供学习交流使用，欢迎大家评论留言或者私信交流，请勿用于任何商业行为。其中，各章实验数据结果也在资源链接中，大家可以自行实验对照，欢迎大家批评指正！文件中每个算法均由作者进行初步测试，部分算法可直接用于工程开发作为模块参考。

各章算法源代码（C++）

Dropdrag

关注

1
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
分类搜索方法与导数方法——《数值计算方法》

数值优化方法的主要内容是利用高效率、高精度的算法求解多元函数极值点的算法。常用的数值优化算法主要包括有分类搜索法、二次逼近方法、内德-米德方法、鲍威尔方法、最速下降法和牛顿方法。首先介绍不同数值优化算法的原理核心思想以及推导过程，分析不同算法的适用场景以及限制条件。再具体分析各个算法模块的基本流程以及核心算法模块架构，最后对每个算法模块进行一到两个实例的数值实验，比较不同算法的精度，分析并提出优缺点以及改进方向。
复制链接

扫一扫