1 scipy.optimize简介

该scipy.optimize包提供几种常用的优化算法。

该模块包含：

1、使用多种算法（例如BFGS，Nelder-Mead单形，牛顿共轭梯度，COBYLA或SLSQP）对多元标量函数进行无约束和无约束的最小化（最小化）

2、全局（强力）优化例程（例如，盆地跳动，differential_evolution）

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

4、标量单变量函数最小化器（minimum_scalar）和根查找器（牛顿）

5、使用多种算法（例如，混合鲍威尔，莱文贝格-马夸特或大型方法，例如牛顿-克里洛夫）的多元方程组求解器（root）。

详见：

#  scipy.optimize import minimizescipy.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 npfrom 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算法

拟牛顿法的核心思想是构造目标函数二阶导数矩阵黑塞矩阵的逆的近似矩阵，避免了解线性方程组求逆的大量计算，更加高效。

介绍：

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, derres = minimize(rosen_and_der, x0, method='BFGS', jac=True,               options={'disp': True})

2.3 牛顿法:Newton-CG

利用黑塞矩阵和梯度来优化。

介绍：

构造目标函数的近似二次型（泰勒展开）：

利用黑塞矩阵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 Hres = 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 Hpres = 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 exampleres = minimize(rosen, x0, method='trust-ncg',               jac=rosen_der, hess=rosen_hess,               options={'gtol': 1e-8, 'disp': True})# Hessian product exampleres = 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 Boundsbounds = Bounds([0, -0.5], [1.0, 2.0])

定义线性约束

from scipy.optimize import LinearConstraintlinear_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 NonlinearConstraintnonlinear_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 varyOptimization terminated successfully.    (Exit mode 0)            Current function value: 0.342717574857755            Iterations: 5            Function evaluations: 6            Gradient evaluations: 5print(res.x)# [0.41494475 0.1701105 ]

3.3 最小二乘:least_squares

这里的f表示损失函数。

假设其为：

其导数：

约束：

from scipy.optimize import least_squaresdef 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) - ydef 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 Ju = 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 pltu_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_scalarf = lambda x: (x - 2) * (x + 1)**2res = minimize_scalar(f, method='brent')print(res.x)# 1.0

5 有界最小化

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

6 自定义最小化器

from scipy.optimize import OptimizeResultdef 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 rootdef 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, dfsol = root(func2, [1, 1], jac=True, method='lm')sol.x# array([ 6.50409711,  0.90841421])

8 线性规划

问题：

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

定义 x：

则目标函数的系数为：

考虑两个不等式约束条件。

第一个是“小于”不等式，所以它已经是linprog所接受的形式。第二个是“大于”不等式，所以需要在两边同时乘以-1，将其转化为“小于”不等式。

准换成矩阵形式：

考虑两个等式：

矩阵形式：

from scipy.optimize import linprogc = 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 Nonex3_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 varymessage: 'Optimization terminated successfully.'    nit: 4  # may vary  slack: array([ 6.52747190e-10, -2.26730279e-09])  # may vary status: 0success: True      x: array([ 9.41025641,  5.17948718, -0.25641026,  1.64102564])  # may varyx = 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_assignmentrow_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

参考：