Python機器學習（六十二）SciPy 優化

優化是指在某些約束條件下，求解目標函數最優解的過程。機器學習、人工智能中的絕大部分問題都會涉及到求解優化問題。

SciPy的optimize模塊提供了許多常用的數值優化算法，一些經典的優化算法包括線性回歸、函數極值和根的求解以及確定兩函數交點的坐標等。

導入scipy.optimize模塊，如下所示：

from scipy import optimize

標量函數極值求解

fmin_bfgs方法

函數f(x)是一個拋物線，求它的極小值：

 

 

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)g(x)=x2+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

解方程組

如果要對如下方程組求解：

 

 func可以定義為：

def func(x):
    u1,u2,u3 = x
    return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]

示例

解方程組：

 

 

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 樣本數據

我們將使用下面的函數來演示曲線擬合：

 

 示例

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()

生成

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()

輸出