優化器scipy.optimize參考指南

1 scipy.optimize簡介

該scipy.optimize包提供幾種常用的優化算法。

該子產品包含：

1、使用多種算法（例如BFGS，Nelder-Mead單形，牛頓共轭梯度，COBYLA或SLSQP）對多元标量函數進行無限制和無限制的最小化（最小化）

2、全局（強力）優化例程（例如，盆地跳動，differential_evolution）

3、最小二乘最小化（least_squares）和曲線拟合（curve_fit）算法

4、标量單變量函數最小化器（minimum_scalar）和根查找器（牛頓）

5、使用多種算法（例如，混合鮑威爾，萊文貝格-馬誇特或大型方法，例如牛頓-克裡洛夫）的多元方程組求解器（root）。

詳見：

https://docs.scipy.org/doc/scipy/reference/optimize.html

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

from scipy.optimize import minimize

scipy.optimize.minimize(
        fun,  #可調用的目标函數。
                       x0,  #ndarray,初值。（n,）
                       args=(), #額外的參數傳遞給目标函數及其導數
                       method=None, #類型的解算器。應該是其中之一:
         #‘Nelder-Mead’、‘Powell’
         #‘CG’、‘BFGS’
         #‘Newton-CG’、‘L-BFGS-B’ 
         #‘TNC’、‘COBYLA’
         #‘SLSQP’、‘dogleg’ 
         #‘trust-ncg’ 

                       jac=None, #目标函數的雅可比矩陣(梯度向量)。
                                 #僅适用于CG, BFGS, Newton-CG, 
                                 #L-BFGS-B, TNC, SLSQP, dogleg，
                                 #trust-ncg。如果jac是一個布爾值，
                                 #且為True，則假定fun将随目标函數傳回
                                 #梯度。如果為False，則用數值方法估計梯
                                 #度。Jac也可以是傳回目标梯度的可調用對
                                 #象。在這種情況下，它必須接受與樂趣相同
                                 #的論點。
                       hess=None, 
                       hessp=None,#目标函數的Hessian(二階導數矩陣)或
                                  #目标函數的Hessian乘以任意向量p。
                                  #僅适用于Newton-CG, dogleg,
                                  #trust-ncg。隻需要給出一個hessp或
                                  #hess。如果提供了hess，則将忽略
                                  #hessp。如果不提供hess和hessp，則用
                                  #jac上的有限差分來近似Hessian積。
                                  #hessp必須計算Hessian乘以任意向量。
               
                       bounds=None, #變量的邊界(僅适用于L-BFGS-B, 
                                    #TNC和SLSQP)。(min, max)
                                    #對x中的每個元素，定義該參數的
                                    #邊界。當在min或max方向上沒有邊界
                                    #時，使用None表示其中之一。
                       constraints=(), #限制定義
                              #(僅适用于COBYLA和SLSQP)
                              # 類型有: ‘eq’ for equality, ‘ineq’ for inequality
                       tol=None, #終止的邊界。
                       callback=None, 
                       options=None)

傳回值: res : OptimizeResult
#以OptimizeResult對象表示的優化結果。重要的屬性有:x是解決方案數組，
#success是一個布爾标志，訓示優化器是否成功退出，以及描述終止原因的消息。
           

2 無限制最小化多元标量函數

2.1 單純形法:Nelder-Mead

函數Rosenbrock ：

x=1時，取最小值。

def rosen(x):
    """The Rosenbrock function"""
    return sum(100.0 * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0)
           

求解：

import numpy as np
from scipy.optimize import minimize

# 初始疊代點
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
 
# 最小化優化器，方法：Nelder-Mead（單純形法）
res = minimize(rosen, x0, method='nelder-mead',
               options={'xatol': 1e-8, 'disp': True})
 
print(res.x)
 
#
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 339
         Function evaluations: 571
[1. 1. 1. 1. 1.]
           

求帶有參數的 Rosenbrock 函數：

def rosen_with_args(x, a, b):
    """The Rosenbrock function with additional arguments"""
    return sum(a*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0) + b
    
res = minimize(rosen_with_args, x0, method='nelder-mead',
               args=(0.5, 1.), options={'xatol': 1e-8, 'disp': True})
           

2.2 拟牛頓法:BFGS算法

拟牛頓法的核心思想是構造目标函數二階導數矩陣黑塞矩陣的逆的近似矩陣，避免了解線性方程組求逆的大量計算，更加高效。

介紹：

https://blog.csdn.net/jiang425776024/article/details/87602847

Rosenbrock導數：

def rosen_der(x):
    # rosen函數的雅可比矩陣
    xm = x[1:-1]
    xm_m1 = x[:-2]
    xm_p1 = x[2:]
    der = np.zeros_like(x)
    der[1:-1] = 200 * (xm - xm_m1 ** 2) - 400 * (xm_p1 - xm ** 2) * xm - 2 * (1 - xm)
    der[0] = -400 * x[0] * (x[1] - x[0] ** 2) - 2 * (1 - x[0])
    der[-1] = 200 * (x[-1] - x[-2] ** 2)
    return der
           

求解：

# 初始疊代點
res = minimize(rosen, x0, method='BFGS', jac=rosen_der,
               options={'disp': True}) 
print(res.x)
           

提供梯度資訊的另一種方法是編寫一個傳回目标和梯度的函數:這可以通過設定jac=True來表示。在這種情況下，要優化的Python函數必須傳回一個元組，其第一個值是目标，第二個值表示梯度。

def rosen_and_der(x):
    objective = sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
    xm = x[1:-1]
    xm_m1 = x[:-2]
    xm_p1 = x[2:]
    der = np.zeros_like(x)
    der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm)
    der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0])
    der[-1] = 200*(x[-1]-x[-2]**2)
    return objective, der

res = minimize(rosen_and_der, x0, method='BFGS', jac=True,
               options={'disp': True})
           

2.3 牛頓法:Newton-CG

利用黑塞矩陣和梯度來優化。

介紹：

https://blog.csdn.net/jiang425776024/article/details/87602847

構造目标函數的近似二次型（泰勒展開）：

利用黑塞矩陣H和梯度做疊代：

黑塞矩陣：

def rosen_hess(x):
    x = np.asarray(x)
    H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
    diagonal = np.zeros_like(x)
    diagonal[0] = 1200*x[0]**2-400*x[1]+2
    diagonal[-1] = 200
    diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
    H = H + np.diag(diagonal)
    return H

res = minimize(rosen, x0, method='Newton-CG',
               jac=rosen_der, hess=rosen_hess,
               options={'xtol': 1e-8, 'disp': True})
res.x
# array([1.,  1.,  1.,  1.,  1.])
           

對于較大的最小化問題，存儲整個Hessian矩陣會消耗大量的時間和記憶體。可将矩陣寫成目标函數的形式：

def rosen_hess_p(x, p):
    x = np.asarray(x)
    Hp = np.zeros_like(x)
    Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1]
    Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1] \
               -400*x[1:-1]*p[2:]
    Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1]
    return Hp

res = minimize(rosen, x0, method='Newton-CG',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'xtol': 1e-8, 'disp': True})

res.x
# array([1., 1., 1., 1., 1.])
           

當Hessian條件不佳時，Newton-CG算法可能是低效的，因為在這些情況下，該方法提供的搜尋方向品質很差。trust-ncg方法可以更有效地處理這種有問題的情況，下面将對此進行描述。

2.4 共轭梯度算法:trust-krylov

與trust-ncg方法類似，trust-krylov方法是一種适用于大規模問題的方法，因為它隻使用hessian作為線性算子，通過矩陣-向量乘積。它比trust-ncg方法更準确地解決了二次子問題。

Newton-CG方法是一種直線搜尋方法：它找到一個搜尋方向，使函數的二次逼近最小化，然後使用直線搜尋算法找到該方向(接近)的最佳步長。另一種方法是，首先固定步長限制，然後通過求解以下二次問題在給定信任半徑内找到最優步長：

根據二次模型與實函數的一緻程度，更新解，調整信任半徑。這類方法稱為信任域方法(trust-region methods)。trust-ncg算法是一種利用共轭梯度算法求解信任域子問題的信任域方法。

# Full Hessian example
res = minimize(rosen, x0, method='trust-ncg',
               jac=rosen_der, hess=rosen_hess,
               options={'gtol': 1e-8, 'disp': True})

# Hessian product example
res = minimize(rosen, x0, method='trust-ncg',
               jac=rosen_der, hessp=rosen_hess_p,
               options={'gtol': 1e-8, 'disp': True})
           

3 限制最小化多元标量函數

最小化函數提供了限制最小化算法，即'trust-constr'， 'SLSQP'和'COBYLA'。它們要求使用稍微不同的結構來定義限制。

方法'trust-constr'要求将限制定義為線性限制和非線性限制對象的序列。

另一方面，方法'SLSQP'和'COBYLA'要求将限制定義為一個字典序列，包含鍵type、fun和jac。

3.1 信任域:trust-constr

信任域限制方法處理的限制最小化問題形式為:

• 定義限制的邊界：

from scipy.optimize import Bounds
bounds = Bounds([0, -0.5], [1.0, 2.0])
           

• 定義線性限制

from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])
           

• 定義非線性限制

def cons_f(x):
    return [x[0]**2 + x[1], x[0]**2 - x[1]]
def cons_J(x):
    return [[2*x[0], 1], [2*x[0], -1]]
def cons_H(x, v):
    return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]])
from scipy.optimize import NonlinearConstraint
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)
           

x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess,
               constraints=[linear_constraint, nonlinear_constraint],
               options={'verbose': 1}, bounds=bounds)

# may vary
# `gtol` termination condition is satisfied.
# Number of iterations: 12, function evaluations: 8, CG iterations: 7, 
# optimality: 2.99e-09, constraint violation: 1.11e-16, execution time: 
# 0.016 s.

print(res.x)

# [0.41494531 0.17010937]
           

3.2 順序最小二乘規劃算法:SLSQP

限制形式：

ineq_cons = {'type': 'ineq',
             'fun' : lambda x: np.array([1 - x[0] - 2*x[1],
                                         1 - x[0]**2 - x[1],
                                         1 - x[0]**2 + x[1]]),
             'jac' : lambda x: np.array([[-1.0, -2.0],
                                         [-2*x[0], -1.0],
                                         [-2*x[0], 1.0]])}
eq_cons = {'type': 'eq',
           'fun' : lambda x: np.array([2*x[0] + x[1] - 1]),
           'jac' : lambda x: np.array([2.0, 1.0])}
           

注意這裡的不等式是≥0，與pymoo(≤0)相反。

x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='SLSQP', jac=rosen_der,
               constraints=[eq_cons, ineq_cons], options={'ftol': 1e-9, 'disp': True},
               bounds=bounds)
# may vary
Optimization terminated successfully.    (Exit mode 0)
            Current function value: 0.342717574857755
            Iterations: 5
            Function evaluations: 6
            Gradient evaluations: 5
print(res.x)
# [0.41494475 0.1701105 ]
           

3.3 最小二乘:least_squares

這裡的f表示損失函數。

假設其為：

其導數：

限制：

from scipy.optimize import least_squares

def model(x, u):
    return x[0] * (u ** 2 + x[1] * u) / (u ** 2 + x[2] * u + x[3])

def fun(x, u, y):
    return model(x, u) - y

def jac(x, u, y):
    J = np.empty((u.size, x.size))
    den = u ** 2 + x[2] * u + x[3]
    num = u ** 2 + x[1] * u
    J[:, 0] = num / den
    J[:, 1] = x[0] * u / den
    J[:, 2] = -x[0] * num * u / den ** 2
    J[:, 3] = -x[0] * num / den ** 2
    return J

u = np.array([4.0, 2.0, 1.0, 5.0e-1, 2.5e-1, 1.67e-1, 1.25e-1, 1.0e-1,
              8.33e-2, 7.14e-2, 6.25e-2])
y = np.array([1.957e-1, 1.947e-1, 1.735e-1, 1.6e-1, 8.44e-2, 6.27e-2,
              4.56e-2, 3.42e-2, 3.23e-2, 2.35e-2, 2.46e-2])
x0 = np.array([2.5, 3.9, 4.15, 3.9])
res = least_squares(fun, x0, jac=jac, bounds=(0, 100), args=(u, y), verbose=1)
# may vary
`ftol` termination condition is satisfied.
Function evaluations 130, initial cost 4.4383e+00, final cost 1.5375e-04, first-order optimality 4.92e-08.

res.x
# array([ 0.19280596,  0.19130423,  0.12306063,  0.13607247])

import matplotlib.pyplot as plt
u_test = np.linspace(0, 5)
y_test = model(res.x, u_test)
plt.plot(u, y, 'o', markersize=4, label='data')
plt.plot(u_test, y_test, label='fitted model')
plt.xlabel("u")
plt.ylabel("y")
plt.legend(loc='lower right')
plt.show()           

4 單變量函數最小化器

from scipy.optimize import minimize_scalar
f = lambda x: (x - 2) * (x + 1)**2
res = minimize_scalar(f, method='brent')
print(res.x)
# 1.0
           

5 有界最小化

from scipy.special import j1
res = minimize_scalar(j1, bounds=(4, 7), method='bounded')
res.x
# 5.33144184241
           

6 自定義最小化器

from scipy.optimize import OptimizeResult
def custmin(fun, x0, args=(), maxfev=None, stepsize=0.1,
        maxiter=100, callback=None, **options):
    bestx = x0
    besty = fun(x0)
    funcalls = 1
    niter = 0
    improved = True
    stop = False

    while improved and not stop and niter < maxiter:
        improved = False
        niter += 1
        for dim in range(np.size(x0)):
            for s in [bestx[dim] - stepsize, bestx[dim] + stepsize]:
                testx = np.copy(bestx)
                testx[dim] = s
                testy = fun(testx, *args)
                funcalls += 1
                if testy < besty:
                    besty = testy
                    bestx = testx
                    improved = True
            if callback is not None:
                callback(bestx)
            if maxfev is not None and funcalls >= maxfev:
                stop = True
                break

    return OptimizeResult(fun=besty, x=bestx, nit=niter,
                          nfev=funcalls, success=(niter > 1))

x0 = [1.35, 0.9, 0.8, 1.1, 1.2]
res = minimize(rosen, x0, method=custmin, options=dict(stepsize=0.05))

res.x
# array([1., 1., 1., 1., 1.])
           

7 找根

• 單變量

from scipy.optimize import root

def func(x):
    return x + 2 * np.cos(x)

sol = root(func, 0.3)

sol.x
# array([-1.02986653])

sol.fun
# array([ -6.66133815e-16])
           

• 多變量

def func2(x):
    f = [x[0] * np.cos(x[1]) - 4,
         x[1]*x[0] - x[1] - 5]
    df = np.array([[np.cos(x[1]), -x[0] * np.sin(x[1])],
                   [x[1], x[0] - 1]])
    return f, df

sol = root(func2, [1, 1], jac=True, method='lm')

sol.x
# array([ 6.50409711,  0.90841421])
           

8 線性規劃

問題：

上述是最大小，需要轉換成最小化，才可使用求解器linprog。

定義 x：

則目标函數的系數為：

考慮兩個不等式限制條件。

第一個是“小于”不等式，是以它已經是linprog所接受的形式。第二個是“大于”不等式，是以需要在兩邊同時乘以-1，将其轉化為“小于”不等式。

準換成矩陣形式：

考慮兩個等式：

矩陣形式：

from scipy.optimize import linprog

c = np.array([-29.0, -45.0, 0.0, 0.0])

A_ub = np.array([[1.0, -1.0, -3.0, 0.0],
                [-2.0, 3.0, 7.0, -3.0]])
b_ub = np.array([5.0, -10.0])

A_eq = np.array([[2.0, 8.0, 1.0, 0.0],
                [4.0, 4.0, 0.0, 1.0]])
b_eq = np.array([60.0, 60.0])

x0_bounds = (0, None)
x1_bounds = (0, 5.0)
x2_bounds = (-np.inf, 0.5)  # +/- np.inf can be used instead of None
x3_bounds = (-3.0, None)

bounds = [x0_bounds, x1_bounds, x2_bounds, x3_bounds]

result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds)

print(result)

     con: array([15.5361242 , 16.61288005])  # may vary
     fun: -370.2321976308326  # may vary
 message: 'The algorithm terminated successfully and determined that the problem is infeasible.'
     nit: 6  # may vary
   slack: array([ 0.79314989, -1.76308532])  # may vary
  status: 2
 success: False
       x: array([ 6.60059391,  3.97366609, -0.52664076,  1.09007993])  # may vary
           

結果表明問題是不可行的，這意味着不存在滿足所有限制條件的解向量。

這并不一定意味着做錯。有些問題确實是行不通的。

假設限制太緊了，可以放松。調整代碼x1_bounds =(0,6)：

x1_bounds = (0, 6)
bounds = [x0_bounds, x1_bounds, x2_bounds, x3_bounds]

result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds)

print(result)

    con: array([9.78840831e-09, 1.04662945e-08])  # may vary
    fun: -505.97435889013434  # may vary
message: 'Optimization terminated successfully.'
    nit: 4  # may vary
  slack: array([ 6.52747190e-10, -2.26730279e-09])  # may vary
 status: 0
success: True
      x: array([ 9.41025641,  5.17948718, -0.25641026,  1.64102564])  # may vary

x = np.array(result.x)

print(c @ x)
# -505.97435889013434  # may vary
           

9 指派問題

考慮到為遊泳混合泳接力隊挑選學生的問題。我們有一個表格，列出了五個學生每種泳姿的時間：

解決方案：定義一個成本矩陣C，每行隻有一列有值，矩陣和為最終成本

定義目标函數：

X是一個二值矩陣，當行i被指定到列j時：。

cost = np.array([[43.5, 45.5, 43.4, 46.5, 46.3],
                 [47.1, 42.1, 39.1, 44.1, 47.8],
                 [48.4, 49.6, 42.1, 44.5, 50.4],
                 [38.2, 36.8, 43.2, 41.2, 37.2]])

from scipy.optimize import linear_sum_assignment
row_ind, col_ind = linear_sum_assignment(cost) # 行索引，列索引

row_ind
# array([0, 1, 2, 3])
col_ind
# array([0, 2, 3, 1])

# 最優配置設定
styles = np.array(["backstroke", "breaststroke", "butterfly", "freestyle"])[row_ind]
students = np.array(["A", "B", "C", "D", "E"])[col_ind]
dict(zip(styles, students))
{'backstroke': 'A', 'breaststroke': 'C', 'butterfly': 'D', 'freestyle': 'B'}

# 最少時間
cost[row_ind, col_ind].sum()
# 163.89999999999998
           

參考：

https://blog.csdn.net/xu624735206/article/details/117320847

https://blog.csdn.net/jiang425776024/article/details/87885969

https://docs.scipy.org/doc/scipy/tutorial/optimize.html