计算数学
数值优化
一维搜索之连续抛物线插值
在黄金分割法中，除了比较函数值大小之外，函数值本身并未在最优求解中用到。 这说明，我们并未充分利用现有的信息来加速迭代收敛。这一节介绍连续抛物线插值法 求解一维搜索问题，这是一种较少浪费函数值的方法。下面阐述连续抛物线插值的基本思想，先在最小值附件选3个点$r,s,t$。计算这3个点的目标函数值$f$，且可以 给出过这3个点的抛物线 \[P\left( x \right) = f\left( r \right) + {d_1}\left( {x - r} \right) + {d_3}\left( {x - r} \right)\left( {x - s} \right)\] 其中 \[{d_1} = \frac{{f\left( s \right) - f\left( r \right)}}{{s - r}}\] \[{d_2} = \frac{{f\left( t \right) - f\left( s \right)}}{{t - s}}\] \[{d_3} = \frac{{{d_2} - {d_1}}}{{t - r}}\] 令\(P\left( x \right)\)可以给出抛物线的极小值公式 \[x = \frac{{r + s}}{2} - \frac{{\left( {f\left( s \right) - f\left( r \right)} \right)\left( {t - r} \right)\left( {t - s} \right)}}{{2\left[ {\left( {s - t} \right)\left( {f\left( t \right) - f\left( s \right)} \right) - \left( {f\left( s \right) - f\left( r \right)} \right)\left( {t - s} \right)} \right]}}\] 此点即为极小值新的近似。有些教材给的算法到此为止，这样的话只能给出较粗的极小值。 在实际应用时，可以重复以上过程，用新的x代替r，s，t中最远的或次最远的点， 反复迭代，即可以给出需求精度下的极小值。
下面给出连续抛物线求一维搜索算法：
- [Input] 最大允许迭代次数、误差容限，初始猜想极小点r，s，t 。
- [Output] 极值点的自变量与函数值。
- [STEP1] 对\(i = 1,2,3, \cdots \)，作如下循环 \[x = \frac{{r + s}}{2} - \frac{{\left( {f\left( s \right) - f\left( r \right)} \right)\left( {t - r} \right)\left( {t - s} \right)}}{{2\left[ {\left( {s - t} \right)\left( {f\left( t \right) - f\left( s \right)} \right) - \left( {f\left( s \right) - f\left( r \right)} \right)\left( {t - s} \right)} \right]}}\] \[t = s\] \[s = r\] \[r = x\] 如果$i$超过与设定的最大允许迭代次数，停止循环，并给出告警信息。 如果已经满足预设的精度需求，则停止循环。
- [STEP2] 输出计算结果，并停止该算法。
package main
import (
"fmt"
."math"
)
func main() {
/*------------------------------------------------------
Author : Song Yezhi
verison : 2021-10-7 18:25
go build -gcflags "-N -l"
------------------------------------------------------*/
r :=0.0
s :=0.7
t :=1.0
_ , _ = paraInterp(r,s,t,myfx)
}
func myfx(x float64) float64 {
/*------------------------------------------------------
Author : Song Yezhi
verison : 2021-10-7 18:26
go build -gcflags "-N -l"
------------------------------------------------------*/
fx := Pow(x,6)-11.0*Pow(x,3)+17.0*x*x - 7.0*x +1.0
return fx
}
func paraInterp(r,s,t float64,funcX func(float64) float64) (x, fx float64) {
/*------------------------------------------------------
Author : Song Yezhi
verison : 2021-10-7 18:28
go build -gcflags "-N -l"
连续抛物型法一维搜索
-----------------------------------------------------
Input Parameters :
r,s,t ------初值
funcX ----- 目标函数
Output Parameters :
x---极小点
fx--极小点函数值
------------------------------------------------------*/
var imax int = 20
var tol float64 = 1e-8
var fs,ft,fr float64
var tmp1,tmp2 float64
fs = funcX(s)
fr = funcX(r)
ft = funcX(t)
for i := 0; i < imax; i++ {
tmp1 = (fs -fr) * (t-r) * (t-s)
tmp2 = 2.0*((s-r)*(ft-fs) -(fs-fr)*(t-s))
x = (r+s)/2.0 - tmp1/tmp2
if Abs(fs-fr) < tol {
break
}
t = s
s = r
r = x
ft = fs
fs = fr
fr = funcX(r)
// *+++single function evaluation +++*
fmt.Printf("i= %4d x= %12.7f f(x)=%12.7f \n", i, x, ft)
}
return x,fx
}