数值算法：无约束优化之一维搜索方法之黄金分割法斐波那契数列法
Posted 宏斌PKUCIS
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了数值算法：无约束优化之一维搜索方法之黄金分割法斐波那契数列法相关的知识，希望对你有一定的参考价值。
目标函数为一元单值函数f：R->R的最小化优化问题，一般不会单独遇到，它通常作为多维优化问题中的一个部分出现，例如梯度下降法中每次最优迭代步长的估计。
一维搜索方法是通过迭代方式求解的，这不同于我们人脑的直接通过解表达式求解方法。迭代算法是从初始搜索点x(0)出发，产生一个迭代序列x(1),x(2),...。在第k=0,1,2，...次迭代中，通过当前迭代点x(k)和目标函数 f 构建下一个迭代点x(k+1)。某些算法可能只需要用到迭代点处的目标函数值，而另一些算法还可能用到目标函数的导数 f\'(x)，甚至是二阶导数 f\'\'(x)。
1、问题描述：
如果一元单值函数 f:R->R 在闭区间[a0,b0]上是单峰的， 求解 f 在该区间上的极小点。
2、计算机求解方法：
以下方法的本质是：区间压缩(截取)法，通过不断缩小极小点所在的区间长度到足够的精度水平来逼近极小点。
（1）黄金分割法（每次区间截取比例为固定值）
第一步：初始区间[a0,b0]，设定截取比例为 r，则第一次挑选两个中间点 a1和 b1，它们满足对称要求：a1-a0=b0-b1= r(b0-a0)，显然 r<1/2，如果 f(a1) < f(b1)，那么极小点应该在区间 [a0,b1] 中; 否则，如 f(a1) >= f(b1)，极小点应该位于区间 [a1,b0] 中。
第二步：经过第一次压缩后，极小点所在区间变为[a0,b1]，或者[a1,b0]。假定为[a0,b1]，下面需要第二次挑选中间点 a2和 b2 。为了充分利用第一次的挑选点信息，减少计算次数，那么第二次挑选的点 b2 可以取第一次的挑选点 a1 ，这样就只需要计算b2（即a1）在新区间上的对称点 a2 和 f(a2) 。为了实现这样的想法，那么必须满足下面的要求: r(b1-a0)= b1 - b2=b1-a1 <=> r[b0-r(b0-a0)-a0]=(1-2r)(b0-a0)<=>r2-3r+1=0<=>1-r=(sqrt(5)-1)/2 = 0.618， 即每次截取后保留比例为0.618 故称此方法为黄金分割法。因此 N 步之后，区间长度为（1-r）N=(0.618)N, 即总压缩比。
挑选 a1和 b1 为中间点
假定第一次截取[b1,b0], 然后在 [a0,b1] 上挑选 a2和 b2，为了减少计算以 a1 取代 b2， 只需计算 a2 。
1 #ifndef _GOLDENSECTION_H_ 2 #define _GOLDENSECTION_H_ 3 4 typedef float (* PtrOneVarFunc)(float x); 5 6 void GoldenSectionMethod(float a, float b, PtrOneVarFunc f, float epsilon); 7 8 #endif
#include <iostream> #include<cmath> #include "goldensection.h" using namespace std; void GoldenSectionMethod(float a, float b, PtrOneVarFunc f, float epsilon) { float sectionRatio = 1.5-sqrt(5)/2; float a0,b0,a1,b1; int times; a0 = a; b0 = b; times = 0; while(abs(a0-b0) > epsilon) { times++; #ifdef _DEBUG cout<<times<<"th iteration:f("<<a0<<")="<<f(a0)<<",f("<<b0<<")="<<f(b0)<<endl; #endif a1 = sectionRatio*(b0-a0)+a0; b1 = b0-sectionRatio*(b0-a0); if(f(a1) < f(b1)) { b0 = b1; } else { a0 = a1; } } cout<<"minimum point x:"<<a0<<","<<b0<<endl; }
1 #include<iostream> 2 #include "goldensection.h" 3 4 float FunctionofOneVariable(float x) 5 { 6 return 7*x*x-5*x+2; 7 } 8 9 int main() 10 { 11 GoldenSectionMethod(-50, 50, FunctionofOneVariable, 0.001); 12 return 0; 13 }
1 1th iteration:f(-50)=17752,f(50)=17252 2 2th iteration:f(-11.8034)=1036.26,f(50)=17252 3 3th iteration:f(-11.8034)=1036.26,f(26.3932)=4746.24 4 4th iteration:f(-11.8034)=1036.26,f(11.8034)=918.225 5 5th iteration:f(-2.7864)=70.2803,f(11.8034)=918.225 6 6th iteration:f(-2.7864)=70.2803,f(6.23059)=242.589 7 7th iteration:f(-2.7864)=70.2803,f(2.78641)=42.4164 8 8th iteration:f(-0.65778)=8.31761,f(2.78641)=42.4164 9 9th iteration:f(-0.65778)=8.31761,f(1.47084)=9.78945 10 10th iteration:f(-0.65778)=8.31761,f(0.657782)=1.73983 11 11th iteration:f(-0.15528)=2.94518,f(0.657782)=1.73983 12 12th iteration:f(0.155282)=1.39238,f(0.657782)=1.73983 13 13th iteration:f(0.155282)=1.39238,f(0.465844)=1.18985 14 14th iteration:f(0.273906)=1.15564,f(0.465844)=1.18985 15 15th iteration:f(0.273906)=1.15564,f(0.39253)=1.11591 16 16th iteration:f(0.319217)=1.11721,f(0.39253)=1.11591 17 17th iteration:f(0.34722)=1.10783,f(0.39253)=1.11591 18 18th iteration:f(0.34722)=1.10783,f(0.375223)=1.10943 19 19th iteration:f(0.34722)=1.10783,f(0.364527)=1.10752 20 20th iteration:f(0.353831)=1.10722,f(0.364527)=1.10752 21 21th iteration:f(0.353831)=1.10722,f(0.360441)=1.10722 22 22th iteration:f(0.356356)=1.10715,f(0.360441)=1.10722 23 23th iteration:f(0.356356)=1.10715,f(0.358881)=1.10716 24 24th iteration:f(0.356356)=1.10715,f(0.357916)=1.10715 25 minimum point x:0.356952,0.357916
（2）斐波那契数列法
在黄金分割法进行区间压缩的过程中，参数 r 始终保持不变。如果在允许参数 r 不断调整，比如 第 k 次迭代使用参数 rk，第k+1次迭代使用参数rk+1,以此类推。同时为了减少计算次数，要求每次迭代中只需要计算一次目标函数值（即在每次挑选中间点时，重复利用上一次挑选的一个中间点信息）。那么可以得出下面的等式：
rk+1(1-rk)=1-2rk <=>rk+1=1 - rk/(1-rk)
显然很多序列{rk}，0= <rk<=1/2满足上式要求，比如黄金分割法中的序列。因此需要寻找最优的一种序列，使得最终得到的总压缩比达到最小，即(1-r1)(1-r2)...(1-rN)最小。经过数学分析，发现利用斐波那契数列构造数列{rk}={1-FN-k+1/FN-k+2}，可以满足上面的要求，其中Fk是斐波那契数列第k个元素，k<=N，N步迭代后使得总压缩比为1/FN+1。
但是可以发现，最后一步迭代rN=1/2, 这就意味着两个中间点aN=bN , 它们位于区间中点上，相互重合，这样将使得区间无法再继续压缩下去。因此需要我们人为调整 rN=1/2 - ε, ε是一个很小的正实数。这样最终压缩比可以表示为（1+2*ε）/FN+1。
算法要求预先设定参数 ε、最终区间长度，然后求出迭代次数N，FN+1≥(b0-a0)(1+2ε)/最终区间长度。注意最后一次迭代使用的rN=1/2 - ε！
1 #ifndef _FIBONACCISEQUENCE_H_ 2 #define _FIBONACCISEQUENCE_H_ 3 4 typedef float (* PtrOneVarFunc)(float x); 5 6 7 void FibonacciSequenceMethod(float a, float b, PtrOneVarFunc fi, float epsilon); 8 9 #endif
1 #include<iostream> 2 #include<cmath> 3 #include<vector> 4 #include<cassert> 5 #include "FibonacciSequence.h" 6 7 using namespace std; 8 9 vector<int> fibonacciseq; 10 11 void FibonacciSequenceMethod(float a, float b, PtrOneVarFunc f, float epsilon) 12 { 13 float a0,b0,a1,b1,r,fb_threshold,ee; 14 int n,k,fbi_1,fbi; 15 a0 = a; 16 b0 = b; 17 n = k = 0; 18 19 ee = 0.0001; 20 r=0; 21 22 fb_threshold = (b0-a0)*(1+2*ee)/epsilon; 23 fbi=1; 24 fbi_1=0; 25 26 fibonacciseq.push_back(fbi); 27 28 while(fbi<fb_threshold) 29 { 30 int temp; 31 temp = fbi + fbi_1; 32 fbi_1 = fbi; 33 fbi = temp; 34 fibonacciseq.push_back(fbi); 35 36 n++; 37 } 38 39 n--; 40 41 42 for( k=1; k<n+1; k++) 43 { 44 45 if(k==n) 46 { 47 r = 1- fibonacciseq[1]*1.0f/fibonacciseq[2] - ee; 48 } 49 else 50 { 51 r = 1 - fibonacciseq[n-k+1]*1.0f/fibonacciseq[n-k+2]; 52 } 53 54 55 a1 = r*(b0-a0) + a0; 56 b1 = b0 - r*(b0-a0); 57 58 if( f(a1) < f(b1) ) 59 { 60 b0 = b1; 61 } 62 else 63 { 64 a0 = a1; 65 } 66 67 #ifdef _DEBUG 68 cout<<k<<"th iteration:a0="<<a0<<",f(a0)="<<f(a0)<<",b0="<<b0<<",f(b0)="<<f(b0)<<",r="<<r<<endl; 69 #endif 70 } 71 72 cout<<"minimum point x:"<<a0<<","<<b0<<endl; 73 }
1 #include<iostream> 2 #include "FibonacciSequence.h" 3 4 5 float FunctionofOneVariable(float x) 6 { 7 return 7*x*x-5*x+2; 8 } 9 10 int main() 11 { 12 FibonacciSequenceMethod(-50, 50, FunctionofOneVariable, 0.001); 13 return 0; 14 }
1th iteration:a0=-11.8034,f(a0)=1036.26,b0=50,f(b0)=17252,r=0.381966 2th iteration:a0=-11.8034,f(a0)=1036.26,b0=26.3932,f(b0)=4746.24,r=0.381966 3th iteration:a0=-11.8034,f(a0)=1036.26,b0=11.8034,f(b0)=918.225,r=0.381966 4th iteration:a0=-2.7864,f(a0)=70.2803,b0=11.8034,f(b0)=918.225,r=0.381966 5th iteration:a0=-2.7864,f(a0)=70.2803,b0=6.23059,f(b0)=242.589,r=0.381966 6th iteration:a0=-2.7864,f(a0)=70.2803,b0=2.78641,f(b0)=42.4164,r=0.381966 7th iteration:a0=-0.65778,f(a0)=8.31762,b0=2.78641,f(b0)=42.4164,r=0.381966 8th iteration:a0=-0.65778,f(a0)=8.31762,b0=1.47084,f(b0)=9.78945,r=0.381966 9th iteration:a0=-0.65778,f(a0)=8.31762,b0=0.657782,f(b0)=1.73983,r=0.381966 10th iteration:a0=-0.15528,f(a0)=2.94518,b0=0.657782,f(b0)=1.73983,r=0.381966 11th iteration:a0=0.155282,f(a0)=1.39238,b0=0.657782,f(b0)=1.73983,r=0.381966 12th iteration:a0=0.155282,f(a0)=1.39238,b0=0.465843,f(b0)=1.18985,r=0.381967 13th iteration:a0=0.273905,f(a0)=1.15564,b0=0.465843,f(b0)=1.18985,r=0.381963 14th iteration:a0=0.273905,f(a0)=1.15564,b0=0.392528,f(b0)=1.11591,r=0.381974 15th iteration:a0=0.319212,f(a0)=1.11721,b0=0.392528,f(b0)=1.11591,r=0.381944 16th iteration:a0=0.347221,f(a0)=1.10783,b0=0.392528,f(b0)=1.11591,r=0.382022 17th iteration:a0=0.347221,f(a0)=1.10783,b0=0.375229,f(b0)=1.10943,r=0.381818 18th iteration:a0=0.347221,f(a0)=1.10783,b0=0.36452,f(b0)=1.10752,r=0.382353 19th iteration:a0=0.353811,f(a0)=1.10722,b0=0.36452,f(b0)=1.10752,r=0.380952 20th iteration:a0=0.353811,f(a0)=1.10722,b0=0.360401,f(b0)=1.10722,r=0.384615 21th iteration:a0=0.356282,f(a0)=1.10715,b0=0.360401,f(b0)=1.10722,r=0.375 22th iteration:a0=0.356282,f(a0)=1.10715,b0=0.358753,f(b0)=1.10716,r=0.4 23th iteration:a0=0.356282,f(a0)=1.10715,b0=0.35793,f(b0)=1.10715,r=0.333333 24th iteration:a0=0.357106,f(a0)=1.10714,b0=0.35793,f(b0)=1.10715,r=0.4999 minimum point x:0.357106,0.35793
以上是关于数值算法：无约束优化之一维搜索方法之黄金分割法斐波那契数列法的主要内容，如果未能解决你的问题，请参考以下文章