章節
優化是指在某些約束條件下，求解目標函數最優解的過程。機器學習、人工智能中的絕大部分問題都會涉及到求解優化問題。
SciPy的optimize模塊提供了許多常用的數值優化算法，一些經典的優化算法包括線性回歸、函數極值和根的求解以及確定兩函數交點的坐標等。
導入scipy.optimize模塊，如下所示：
from scipy import optimize
標量函數極值求解
fmin_bfgs方法
函數f(x)是一個拋物線，求它的極小值：
$$
f(x) = x^2 + 2x + 9
$$
我們先畫出函數曲線：
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
# 定義函數
def f(x):
return x**2 + 2*x + 9
# x取值：-10到10之間，間隔0.1
x = np.arange(-10, 10, 0.1)
# 畫出函數曲線
plt.plot(x, f(x))
# plt.savefig('./opt2-1.png') # 保存要顯示的圖片
plt.show()
計算該函數最小值的有效方法之一是使用帶起點的BFGS算法。該算法從參數給定的起始點計算函數的梯度下降，並輸出梯度為零、二階導數為正的極小值。BFGS算法是由Broyden，Fletcher，Goldfarb，Shanno四個人分別提出的，故稱為BFGS校正。
示例
使用BFGS函數，找出拋物線函數的最小值:
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
# 定義函數
def f(x):
return x**2 + 2*x + 9
# x取值：-10到10之間，間隔0.1
x = np.arange(-10, 10, 0.1)
# 畫出函數曲線
plt.plot(x, f(x))
# 第一個參數是函數名，第二個參數是梯度下降的起點。返回值是函數最小值的x值(ndarray數組)
xopt = optimize.fmin_bfgs(f, 0)
xmin = xopt[0] # x值
ymin = f(xmin) # y值，即函數最小值
print('xmin: ', xmin)
print('ymin: ', ymin)
# 畫出最小值的點, s=20設置點的大小，c='r'設置點的顏色
plt.scatter(xmin, ymin, s=20, c='r')
#plt.savefig('./opt3-1.png') # 保存要顯示的圖片
plt.show()
輸出
Optimization terminated successfully.
Current function value: 8.000000
Iterations: 2
Function evaluations: 9
Gradient evaluations: 3
xmin: -1.00000000944232
ymin: 8.0
fmin_bfgs有個問題，當函數有局部最小值，該算法會因起始點不同，找到這些局部最小而不是全局最小。
讓我們看另一個函數，該函數有多個局部最小值：
$$
g(x) = x^2 + 20sin(x)
$$
畫出該函數的曲線：
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
# 定義函數
def g(x):
return x**2 + 20*np.sin(x)
# x取值：-10到10之間，間隔0.1
x = np.arange(-10, 10, 0.1)
# 畫出函數曲線
plt.plot(x, g(x))
# plt.savefig('./opt4-1.png') # 保存要顯示的圖片
plt.show()
函數曲線
可以看到，該函數有多個底部。如果初始值設置不當，獲得的最小值有可能只是局部最小值。
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
# 定義函數
def g(x):
return x**2 + 20*np.sin(x)
# x取值：-10到10之間，間隔0.1
x = np.arange(-10, 10, 0.1)
# 畫出函數曲線
plt.plot(x, g(x))
# 第一個參數是函數名，第二個參數是梯度下降的起點。返回值是函數最小值的x值(ndarray數組)
# 可以看到5.0附近有個局部最小，把初始值設置為7, 返回的應該是這個局部最小值。
xopt = optimize.fmin_bfgs(g, 7)
xmin = xopt[0] # x值
ymin = g(xmin) # y值，即函數最小值
print('xmin: ', xmin)
print('ymin: ', ymin)
# 畫出最小值的點, s=20設置點的大小，c='r'設置點的顏色
plt.scatter(xmin, ymin, s=20, c='r')
#plt.savefig('./opt5-1.png') # 保存要顯示的圖片
plt.show()
輸出
Optimization terminated successfully.
Current function value: 0.158258
Iterations: 3
Function evaluations: 27
Gradient evaluations: 9
xmin: 4.271095444673479
ymin: 0.15825752683190686
可以看到5.0附近有個局部最小，把初始值設置為7, 返回的是這個局部最小值。
如果把初始值設置為0，應該就能返回真正的全局最小值：
basinhopping 方法
對於這種情況，可以使用scipy.optimize提供的basinhopping()方法。該方法把局部優化方法與起始點隨機抽樣相結合，求出全局最小值，代價是耗費更多計算資源。
調用語法
optimize.basinhopping(func, x0)
- func 目標函數
- x0 初始值
示例
import numpy as np
from scipy import optimize
# 定義函數
def g(x):
return x**2 + 20*np.sin(x)
# 求取最小值，初始值為7
ret = optimize.basinhopping(g, 7)
print(ret)
輸出
fun: 0.15825752683178962
lowest_optimization_result: fun: 0.15825752683178962
hess_inv: array([[0.04975718]])
jac: array([4.76837158e-07])
message: 'Optimization terminated successfully.'
nfev: 12
nit: 2
njev: 4
status: 0
success: True
x: array([4.27109533])
message: ['requested number of basinhopping iterations completed successfully']
minimization_failures: 0
nfev: 1641
nit: 100
njev: 547
x: array([4.27109533])
可以看到全局最小值被正確找出：fun: 0.15825752683178962，x: array([4.27109533])
scipy.optimize.brute函數蠻力法也可用於全局優化，但效率較低。蠻力方法的語法為:
scipy.optimize.brute(f, 0)
fminbound
要求取一定范圍之內的函數最小值，可使用fminbound方法。
調用語法
optimize.fminbound(func, x1, x2)
- func 目標函數
- x1, x2 范圍邊界
示例
import numpy as np
from scipy import optimize
# 定義函數
def g(x):
return x**2 + 20*np.sin(x)
# 求取-10到-5之間的函數最小值。full_output=True表示返回詳細信息。
ret = optimize.fminbound(g, -10, -5, full_output=True)
print(ret)
輸出
(-7.068891380019064, 35.82273589215206, 0, 12)
函數最小值是：35.82273589215206，對應的x值：-7.068891380019064。
函數求解
對於一個函數f(x)，當f(x) = 0，求取x的值，即為函數求解。這種情況，可以使用fsolve()函數。
調用語法
optimize.fsolve(func， x0)
- func 目標函數
- x0 初始值
解單個方程
示例
解單個方程:
import numpy as np
from scipy import optimize
# 定義函數
def g(x):
return x**2 + 20*np.sin(x)
# 函數求解
ret = optimize.fsolve(g, 2)
print(ret)
輸出
[0.]
解出的根值是：0
解方程組
如果要對如下方程組求解：
$$
f1(u1,u2,u3) = 0 \
f2(u1,u2,u3) = 0 \
f3(u1,u2,u3) = 0
$$
func可以定義為：
def func(x):
u1,u2,u3 = x
return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]
示例
解方程組：
$$
4x + 9 = 0 \
3y^2 - sin(yz) = 0 \
yz - 2.5 = 0
$$
from scipy.optimize import fsolve
from math import sin,cos
def f(x):
x0 = float(x[0])
x1 = float(x[1])
x2 = float(x[2])
return [
4*x1 + 9,
3*x0*x0 - sin(x1*x2),
x1*x2 - 2.5
]
result = fsolve(f, [1,1,1])
print (result)
輸出
[-0.44664383 -2.25 -1.11111111]
擬合
curve_fit
假設有一批數據樣本，要創建這些樣本數據的擬合曲線/函數，可以使用Scipy.optimize模塊的curve_fit()函數。
調用形式
optimize.curve_fit(func, x1, y1)
- func 目標函數
- x1, y1 樣本數據
我們將使用下面的函數來演示曲線擬合：
$$
f(x) = 50cos(x) + 2
$$
示例
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
# 函數模型用於生成數據
def g(x, a, b):
return a*np.cos(x) + b
# 產生含噪聲的樣本數據
x_data = np.linspace(-5, 5, 100) # 采樣點
y_data = g(x_data, 50, 2) + 5*np.random.randn(x_data.size) # 加入隨機數作為噪聲
# 使用curve_fit()函數來估計a和b的值
variables, variables_covariance = optimize.curve_fit(g, x_data, y_data)
# 輸出結果
print('\n求出的系數a, b: ')
print(variables)
print('\nvariables_covariance: ')
print(variables_covariance)
輸出
求出的系數a, b:
[49.66367999 2.09557981]
variables_covariance:
[[0.55593391 0.10388677]
[0.10388677 0.26071478]]
variables是給定模型的最優參數，variables_covariance可用於檢查擬合情況，其對角線元素值表示每個參數的方差。可以看到我們正確求出了系數值。
繪制曲線：
import matplotlib.pyplot as plt
y = g(x_data, variables[0], variables[1])
plt.plot(x_data, y_data, 'o', color="green", label = "Samples")
plt.plot(x_data, y, color="red", label = "Fit")
plt.legend(loc = "best")
#plt.savefig('./opt10-1.png') # 保存要顯示的圖片
plt.show()
生成
[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-98kPUfcs-1571731570441)(https://www.qikegu.com/wp-content/uploads/2019/07/opt10-1.png)]
leastsq()函數
最小二乘法是非常經典的數值優化算法，通過最小化誤差的平方和來尋找最符合數據的曲線。
optimize模塊提供了實現最小二乘擬合算法的函數leastsq()，leastsq是least square的簡寫，即最小二乘法。
調用形式
optimize.leastsq(func, x0, args=())
- func 計算誤差的函數
- x0 是計算的初始參數值
- args 是指定func的其他參數
示例
import numpy as np
from scipy import optimize
# 樣本數據
X = np.array([160,165,158,172,159,176,160,162,171]
Y = np.array([58,63,57,65,62,66,58,59,62])
# 偏差函數, 計算以p為參數的直線和原始數據之間的誤差
def residuals(p):
k, b = p
return Y-(k*X+b)
# leastsq()使得residuals()的輸出數組的平方和最小，參數的初始值為[1, 0]
ret = optimize.leastsq(residuals, [1, 10])
k, b = ret[0]
print("k = ", k, "b = ", b)
輸出
k = 0.4211697393502931 b = -8.288302606523974
繪制曲線：
import matplotlib.pyplot as plt
#畫樣本點
plt.figure(figsize=(8, 6)) ##指定圖像比例： 8：6
plt.scatter(X, Y, color="green", label="Samples", linewidth=2)
#畫擬合直線
x = np.linspace(150, 190, 100) ##在150-190直接畫100個連續點
y = k*x + b ##函數式
plt.plot(x,y,color="red", label="Fit",linewidth=2)
plt.legend() #繪制圖例
plt.savefig('./opt11-1.png') # 保存要顯示的圖片
plt.show()
輸出
