摘要
在生產(chǎn)過程、科學實驗以及日常生活中，人們總希望用最少的人力、物力、財力和時間去辦更多的事，獲得最大的效益，，所以最優(yōu)化理論和方法日益受到重視。
無約束最優(yōu)化計算方法是數(shù)值計算領域中十分活躍的研究課題之一，快速的求解無約束最優(yōu)化問題，除了自身的重要性以外，還體現(xiàn)在它也構成一些約束最優(yōu)化問題的子問題。因此，對于無約束最優(yōu)化問題，如何快速高精度的求解一直是優(yōu)化工作者十分關心的事。本文驗證求解一維無約束最優(yōu)化問題的三種線性搜索方法，分別是牛頓法、黃金分割法，二次插值法。
問題
min ?f(x)=3x4 - 4x3-12x2
1. 一維牛頓法
1.1. 概念
牛頓迭代法（Newton’s method）又稱為牛頓-拉夫遜方法（Newton-Raphson method）,它是牛頓在17世紀提出的一種在實數(shù)域和復數(shù)域上近似求解方程的方法,其基本思想是利用目標函數(shù)的二次Taylor展開，并將其極小化。牛頓法使用函數(shù) f(x)的泰勒級數(shù)的前面幾項來尋找方程 "f(x)=0" 的根。牛頓法是求方程根的重要方法之一，其最大優(yōu)點是在方程 "f(x)=0" 的單根附近具有平方收斂，而且該法還可以用來求方程的重根、復根，此時非線性收斂，但是可通過一些方法變成線性收斂。
1.2. 牛頓法的幾何解釋
方程 f(x)=0 的根 x*可解釋為曲線y=f(x)軸的焦點的橫坐標。如下圖：
設 x[k]是根 x*的某個近似值，過曲線 y=f(x)上橫坐標為 x[k]的點P[k]引切線，并將該切線與 x軸的交點 的橫坐標 x_{k+1}_作為 x*的新的近似值。鑒于這種幾何背景，牛頓法亦稱為切線法。
1.3. 原理
將f(xk+1)在x=xk處一階泰勒展開：
公式不打了
令函數(shù)趨于0，求與x軸交點：
1.4. 程序主體流程圖
1.5. 運行效果
分別取-1.5,0.4,1.0,1.6,2.8
1.6. 評價
1.6.1. 優(yōu)點
在凸函數(shù)很初始點選取好的情況下，收斂快。
1.6.2. 缺點
? 計算二階導數(shù)，計算量大
? 要求函數(shù)二階可微
? 收斂性與初始點選取有關
2. 黃金分割法
2.1. 基本思路
黃金分割法適用于[a,b]區(qū)間上的任何單股函數(shù)求極小值問題，對函數(shù)除要求“單峰”外不做其他要求，甚至可以不連續(xù)。因此，這種方法的適應面非常廣。黃金分割法也是建立在區(qū)間消去法原理基礎上的試探方法，即在搜索區(qū)間 ？內(nèi)適當插入兩點 ？，并計算其函數(shù)值。 ？將區(qū)間分成三段，應用函數(shù)的單峰性質(zhì)，通過函數(shù)值大小的比較，刪去其中一段，是搜索區(qū)間得以縮小。然后再在保留下來的區(qū)間上作同樣的處理，如此迭代下去，是搜索區(qū)間無限縮小，從而得到極小點的數(shù)值近似解。
2.2. 基本原理與步驟
一維搜索是解函數(shù)極小值的方法之一，其解法思想為沿某一已知方向求目標函數(shù)的極小值點。一維搜索的解法很多，這里主要采用黃金分割法（0.618法）。如圖所示
黃金分割法是用于一元函數(shù) ？在給定初始區(qū)間 ？內(nèi)搜索極小點？的一種方法。其基本原理是：依照“去劣存優(yōu)”原則、對稱原則、以及等比收縮原則來逐步縮小搜索區(qū)間。
具體步驟是：在區(qū)間 ？內(nèi)取點： ？分為三段。如果 ？,令？；如果 ？，令 <？,如果 ？都大于收斂精度？重新開始。因為 ？為單峰區(qū)間，這樣每次可將搜索區(qū)間縮小0.618倍或0.382倍，處理后的區(qū)間都將包含極小點的區(qū)間縮小，然后在保留下來的區(qū)間上作同樣的處理，如此迭代下去，將使搜索區(qū) ？逐步縮小，滿足預先給定的精度時，即獲得一維優(yōu)化問題的近似最優(yōu)解。
2.3. 算法流程圖
2.4. 運行效果
[-2,0],[1,3]
2.5. 評價
很巧妙地使每次分割點都為上次的黃金分割點，可以復用上次運算的結果，簡化計算。它是優(yōu)化計算中的經(jīng)典算法，以算法簡單、收斂速度均勻、效果較好而著稱，是許多優(yōu)化算法的基礎，但它只適用于一維區(qū)間上的凸函數(shù)，即只在單峰區(qū)間內(nèi)才能進行一維尋優(yōu)，其收斂效率較低。
3. 二次插值法
3.1. 基本思路
在求解一元函數(shù)?的極小點時，在搜索區(qū)間中用低次（通常不超過三次）插值多項式?來近似目標函數(shù)，后求該多項式的極小點（比較容易計算），并以此作為目標函數(shù)?的近似極小點。如果其近似的程度尚未達到所要求的精度時，反復使用此法，逐次擬合，直到滿足給定的精度時為止。
3.2. 設計思路
考慮二次多項式
則
令 ，得 。這意味著我們要求a，b。
今考慮在包含 的極小點 的搜索區(qū)間 中，給定三個點 ， ， ，滿足
< <
> <
利用三點處的函數(shù)值 ， ， 構造二次函數(shù)，并要求插值條件滿足
令 ，i=1,2,3。解上述方程組得
于是，二次函數(shù) 的極小點為
設 。求得 和 以后，如果
≤ ，當 > 時，
或者如果
≤ ，當 < 時。
則我們認為收斂準則滿足。如果 < ,則極小點估計為 ，否則為 。
若終止準則不滿足，則利用 提供的信息，從 ， ， 和 中選出相鄰的三個點，將原來的搜索區(qū)間縮小，然后重復上述過程，直到終止準則滿足為止。
3.3. 算法設計
初始步 給出?，?，?，滿足上述設計步驟。
步1 由上述設計步驟計算?。
步2 比較?和?的大小，如果?>?，則轉步3；否則轉步4。
步3 如果?≤?，則
?，?，?，?，
轉步5；否則?，?，轉步5。
步4 若?，則
?，?，?，?，
轉步5；否則?，?，轉步5。
步5 如果收斂準則滿足，停止迭代；否則轉步1，在新的搜索區(qū)間[?,?]上按公式計算二次插值函數(shù)的極小點?。
3.4. 程序框圖
3.5. 運行結果
-2.5,-1.4,-0.3
-0.3,1.5,3.3.
3.6. 評價
3.6.1. 優(yōu)點
插值法僅需計算函數(shù)值，不涉及導數(shù)、Hesse矩陣等的計算，計算起來相對比較簡單，能夠適用于非光滑和導數(shù)表達式復雜或表達式寫不出等種種情形。
3.6.2. 缺點
當?shù)綌?shù)較多時，計算過程比較復雜，計算量較大，計算起來比較麻煩。當?shù)c離目標函數(shù)的最優(yōu)解較遠時，追求線性搜索的精度反而會降低整個算法的效率。
程序
用符號推導畫圖
picture.py
import matplotlib.pyplot as plt
import numpy as np
import sympy
t = sympy.symbols('t')
y=3*pow(t,4)-4*pow(t,3)-12*pow(t,2)
y_diff1 = sympy.diff(y, t)
y_diff2 = sympy.diff(y, t, 2)
f=sympy.lambdify(t,y,"numpy")
f1=sympy.lambdify(t,y_diff1,"numpy")
f2=sympy.lambdify(t,y_diff2,"numpy")
# print(sympy.solveset(sympy.Eq(y_diff1,0)))
station=sympy.solveset(sympy.Eq(y_diff1,0))
station=list(station)
station=list(map(int,station))
# 去除鞍點
for i,a in enumerate(station):
if f2(np.array(a))==[0]:
station=station.remove(a)
pass
print("所有駐點",station)
station.sort()
region=0.2*(station[-1]-station[0])
t=np.arange(station[0]-region,station[-1]+region,0.01)
fg,axes=plt.subplots(3)
axes[0].plot(t,f(t))
axes[0].set_title("original function")
axes[1].plot(t,f1(t))
axes[1].plot(station,f1(np.array(station)),'ro')
# axes[1].set_title("一階導數(shù)")
axes[2].plot(t,f2(t))
# axes[2].set_title("二階導數(shù)")
for i in range(3):
axes[i].grid()
pass
plt.show()
Gold_ratio.py
import math
import numpy as np
import matplotlib.pyplot as plt
def draw_dot(x,y,tx):
plt.plot([x],[y],'o')
plt.text(x,y,tx,color='blue',fontsize=15)
pass
Epsilon=1e-2#極限
# y=lambda t:t*(t+2)#函數(shù)
y=lambda t:3*pow(t,4)-4*pow(t,3)-12*pow(t,2)
(a,b)=[-2,0]# 確定a,b
draw_dot(a,y(a),'a')
draw_dot(b,y(b),'b')
Belta=(3-math.sqrt(5))/2 # 0.381
t1=a+Belta*(b-a)
y1=y(t1)
t2=a+(1-Belta)*(b-a)
y2=y(t2)
n=0
t = np.arange(a-0.1*(b-a), b+0.2*(b-a), 0.01)
s = y(t)
plt.plot(t, s)
# draw_dot(t2,y(t2)+1*n,'r'+str(n))
# draw_dot(t1,y(t1)+1*n,'l'+str(n))
while (t2-t1)>Epsilon:
n+=1
if y1>y2:
a=t1
t1=t2
t2=a+b-t1
y1=y2
y2=y(t2)
# draw_dot(a,y(a),'a'+str(n))
# draw_dot(t2,y(t2),'t2+'+str(n))
draw_dot(t2,y(t2)+2*n,'r'+str(n))
draw_dot(t1,y(t1)+2*n,'l'+str(n))
pass
else:
b=t2
t2=t1
y2=y1
t1=a+b-t2
y1=y(t1)
# draw_dot(b,y(b),'b'+str(n))
draw_dot(t2,y(t2)+1*n,'r'+str(n))
draw_dot(t1,y(t1)+1*n,'l'+str(n))
#繪制函數(shù)圖像
t_ = (t1 + t2) / 2
y_ = y(t_)
print('t*=%f,y*=%f' % (t_, y_))
def printFunc(x, y):
plt.plot([x], [y], 'ro')
plt.text(x,y,'(%0.3f,%0.3f)'%(x,y), color='red', fontsize=15)
printFunc(t_,y_)
plt.show()
Newton.py
import sympy
import matplotlib.pyplot as plt
import numpy as np
Epsilon = 0.01
t = sympy.symbols('t')
# y=t*(t+2)
y = t ** 3 - 2 * t + 1
f=lambda t:t ** 3 - 2 * t + 1
T0 = [0.5]#初始值
# y = 3 * pow(t, 4) - 4 * pow(t, 3) - 12 * pow(t, 2)
# f = lambda t: 3 * pow(t, 4) - 4 * pow(t, 3) - 12 * pow(t, 2)
# T0 = [-1.5,0.4,1.0,1.6,2.8]
y_diff1 = sympy.diff(y, t)
y_diff2 = sympy.diff(y, t, 2)
k = 1
def printFunc(f, a, b):
t = np.arange(a, b, 0.01)
s = f(t)
plt.plot(t, s)
plt.show()
def Next_t(t0):
global k
delta = (y_diff1 / y_diff2).subs(t, t0)
t1 = (t - delta).subs(t, t0)
# print(t1)
if abs(delta) < Epsilon:
print('第%d次迭代：%f=%f%+f' % (k, t1, t0, -delta))
plt.plot([t1], [f(t1)], marker='o', markersize=8, markerfacecolor="red")
plt.text(t1, f(t1), k, color='red', fontsize=k * 2 + 10)
return t1
else:
print('第%d次迭代：%f=%f%+f' % (k, t1, t0, -delta))
plt.plot([t1], [f(t1)], marker='o', markersize=8, markerfacecolor="green")
plt.text(t1, f(t1), k, color='green', fontsize=k * 2 + 10)
k += 1
return Next_t(t1)
for t0 in T0:
Next_t(t0)
printFunc(f, -3, 3)
# printFunc(f,t0-1,t0+1)
Polynomial_interpolation.py
'''
二次插值法python實現(xiàn)
f(x)=x^4 - 4x^3 - 6x^2 -16x +4極值
區(qū)間[-1,6] e=0.05
'''
import numpy as np
import matplotlib.pyplot as plt
from time import sleep
from threading import Thread
'''
函數(shù)表達式
'''
def f(x):
# return 1.0 * (pow(x, 4) - 4 * pow(x, 3) - 6 * pow(x, 2) - 16 * x + 4)
return 3*pow(x,4)-4*pow(x,3)-12*pow(x,2)
# 定義變量們
X2,Y=list(),list()
k=0
a = -2.5# 左點
b = -0.3# 右點
e = 0.001# 精度
'''
繪制函數(shù)圖像
'''
def close(time=1):
sleep(time)
# plt.savefig('./img/'+name)
plt.close()
pass
def printFunc():
t = np.arange(a, b, 0.01)
s = f(t)
plt.plot(t, s)
def update_point(x,y):
global k
printFunc()
plt.plot(x,y,'ro')
plt.text(x[-1], y[-1], k, color='red', fontsize=k + 10)
# else:
# plt.plot([x], [y], 'ro')
# plt.text(x, y, k, color='red', fontsize=k + 10)
thread1 = Thread(target=close, args=())
thread1.start()
# print('打開')
plt.show()
# print("close")
def final_fun(x,y):
global k
printFunc()
plt.plot(x, y, 'ro')
for i in range(1,k+1):
plt.text(x[i-1], y[i-1], i, color='red', fontsize=i + 10)
# thread1 = Thread(target=close, args=())
# thread1.start()
plt.show()
'''
e為精度
'''
def search(f, x1, x2, x3):
global k
k+=1
if f(x2) > f(x1) or f(x2) > f(x3):
print("不滿足兩頭大中間小的性質(zhì)")
return 0
# 系數(shù)矩陣
A = [[pow(x1, 2), x1, 1], [pow(x2, 2), x2, 1], [pow(x3, 2), x3, 1]]
b = [f(x1), f(x2), f(x3)]
X = np.linalg.solve(A, b)
a0,a1,_ = X
x = - a1 / (2 * a0)
# 達到精度退出
if abs(x - x2) < e:
if f(x) < f(x2):
y = f(x)
print('最后的x:',x)
X2.append(x)
Y.append(y)
final_fun(X2,Y)
return (X2, Y)
else:
y = f(x2)
print("最后的x2",x2)
X2.append(x2)
Y.append(y)
final_fun(X2, Y)
return (X2, Y)
arr = [x1, x2, x3, x]
arr.sort()
# 在x2和新算出的x中找最小值
if f(x2) > f(x):
index = arr.index(x)
x2 = x
else:
index = arr.index(x2)
x1 = arr[index - 1]
x3 = arr[index + 1]
X2.append(x2)
Y.append(f(x2))
print('運行中的第%d次：%f' % (k, x2))
update_point(X2,Y)
return search(f, x1, x2, x3)
def regre(f, a, b):
x1 = a
x3 = b
x2 = (a + b) / 2.0
search(f, x1, x2, x3)
regre(f, a, b)
