SciPy 最优化之最小化 | 隔叶黄莺 Yanbin Blog

2024-12-05 | 阅读(14)

SciPy 是一个开源的算法库和数学工具包，可以处理最优化、线性代数、积分、插值、拟合、特殊函数、快速傅里叶变换、信号处理、图像处理、常微分方程求解器等。它依赖于 NumPy, Pandas 也依赖了 NumPy。本文重点是体验它怎么处理最优化的问题。很多情形下通过 SciPy 的 optimize.minimize 方法寻求目标函数最小值的过程得到最优化的输入与输出。比如寻找二次元函数的根，求解线性/动态规则，金融行业的计算出最优投资组合的资产分配等。为什么 SciPy 没有 maximize 方法呢，因为没有必要，想要找到最大化的值，只要把目标函数的值取反，或者是模或绝对值的最小值。看到 minimize 方法名更让人觉得目标函数会有一个收敛值。

虽然 SciPy 对特定的问题有更直白的函数，如求根有 optimize.root, 线性规则 optimize.linprog(现不建议使用)，但各种优化基本都可以回归到 minimize 方法调用。minimize 方法的原型是

def minimize(fun, x0, args=(), method=None, jac=None, hess=None,
             hessp=None, bounds=None, constraints=(), tol=None,
             callback=None, options=None):

def minimize(fun, x0, args=(), method=None, jac=None, hess=None,

hessp=None, bounds=None, constraints=(), tol=None,

callback=None, options=None):

除了必须的目标函数和初始值，还有更多参数，像常用的约束(contraints) - 满足某些特定条件的最优化, 线程或非线性约束等; 求解方法(method) - Powell, Newton-CG 等

下面用 optimize.minimize 来求解一些问题

求一元二次函数的最小值(根)

函数定义为 $f(x) = 4x^2 + 2x + 10$ , 画出它的函数曲线如下

根(顶点) 的 x 值是 $x = -\frac{b}{2a} = -\frac{2}{8} = -0.25$ , 相应的 y 值是 9.75

现在来看用 minimize 怎么求解的，代码如下

import scipy.optimize as opt
import numpy as np

def objective(x):
    return 4*(x**2)+2*x+10

x0 =  np.array([-1])
result = opt.minimize(objective, x0)

print(result.x)           # 目标函数最小值时输入 x
print(result.fun)         # 目标函数最小值
print(result)             # 迭代结果

import scipy.optimize as opt

import numpy as np

def objective(x):

return 4*(x**2)+2*x+10

x0 = np.array([-1])

result = opt.minimize(objective, x0)

print(result.x) # 目标函数最小值时输入 x

print(result.fun) # 目标函数最小值

print(result) # 迭代结果

x0 是猜测的初始值，minimize 通过迭代算出目标函数最小值的。

上面代码输出为

[-0.25]
9.75
  message: Optimization terminated successfully.
  success: True
   status: 0
      fun: 9.75
        x: [-2.500e-01]
      nit: 2
      jac: [ 0.000e+00]
 hess_inv: [[ 1.250e-01]]
     nfev: 6
     njev: 3

[-0.25]

9.75

message: Optimization terminated successfully.

success: True

status: 0

fun: 9.75

x: [-2.500e-01]

nit: 2

jac: [ 0.000e+00]

hess_inv: [[ 1.250e-01]]

nfev: 6

njev: 3

函数最小值为 9.75, x 为 -0.25，通过两次(nit) 迭代算出。

不同的初始值对结果也会影响

比如 x0 = np.array[10] 时， x=-0.2500001, f(x)=9.750000000000037

-1 和 10 是 -0.25 两边的值，所以 minimize 会向两边进行迭代，不同初始值会影响到迭代的次数。

我们可以输出它的迭代过程，用 callback 函数

def callback(xk):
    print(f"Iteration: {callback.iter}, x = {xk}, y = {objective(xk)}")
    callback.iter += 1

callback.iter = 0

x0 =  np.array([100])
result = opt.minimize(objective, x0, callback=callback)

def callback(xk):

print(f"Iteration: {callback.iter}, x = {xk}, y = {objective(xk)}")

callback.iter += 1

callback.iter = 0

x0 = np.array([100])

result = opt.minimize(objective, x0, callback=callback)

输出为

Iteration: 0, x = [78.79], y = [24999.0364]
Iteration: 1, x = [30.19728879], y = [3717.89957945]
Iteration: 2, x = [-0.24999419], y = [9.75]
Iteration: 3, x = [-0.25], y = [9.75]

Iteration: 0, x = [78.79], y = [24999.0364]

Iteration: 1, x = [30.19728879], y = [3717.89957945]

Iteration: 2, x = [-0.24999419], y = [9.75]

Iteration: 3, x = [-0.25], y = [9.75]

很神奇，不知道它内部是怎么算定出 x=-0.25 时就是最小值的

如果是一个下开口的二次函数 $f(x) = -4x^2 + 2x + 10$ , 也试图求其最小值会如何呢？

def objective(x):
    return -4*(x**2)+2*x+10

x0 =  np.array([10])
result = opt.minimize(objective, x0)
print(result)

def objective(x):

return -4*(x**2)+2*x+10

x0 = np.array([10])

result = opt.minimize(objective, x0)

print(result)

输出结果

  message: Desired error not necessarily achieved due to precision loss.
  success: False
   status: 2
      fun: -4359650.2304
        x: [ 1.044e+03]
      nit: 1
      jac: [-8.352e+03]
 hess_inv: [[-1.250e-01]]
     nfev: 236
     njev: 112

message: Desired error not necessarily achieved due to precision loss.

success: False

status: 2

fun: -4359650.2304

x: [ 1.044e+03]

nit: 1

jac: [-8.352e+03]

hess_inv: [[-1.250e-01]]

nfev: 236

njev: 112

迭代一次就马上止住，给出答案为没有最小值。但只要改造一下目标函数, 即取反

def objective(x):
    return -(-4*(x**2)+2*x+10)

1 2	def objective(x): return -(-4(x2)+2x+10)

这样就能算出它 x= 0.25 时最大值为 10.25，因为向下开口图形的顶点为最大值，取反则为最小值。

我们知道直线没有最小最大值，用 SciPy 试下求函数 $f(x) = 4x + 10$ 的最小值

def objective(x):
    return 4 * x + 10

x0 = np.array([10])
result = opt.minimize(objective, x0)

def objective(x):

return 4 * x + 10

x0 = np.array([10])

result = opt.minimize(objective, x0)

果然，结局是

message: Desired error not necessarily achieved due to precision loss.
success: False

求立体空间中的最小值

如果是一个二元的函数 $f(x, y) = 4x^2 + 2y^2 + 10$ , 我们先画出它的图形来

现在是看图索骥了，就是想要知道 x, y 分别为多少时函数最小。把这个图形看成是一个地势高低图的话，就是在寻找畦地

def objective(x_arr):
    x, y = x_arr
    return 4*(x**2) + 2*(y**2) + 10

x0 = np.array([-10, -10])
result = opt.minimize(objective, x0)

np.set_printoptions(suppress=True)  # 数值不以科学计数法输出
print(result.x)
print(result.fun)

def objective(x_arr):

x, y = x_arr

return 4*(x**2) + 2*(y**2) + 10

x0 = np.array([-10, -10])

result = opt.minimize(objective, x0)

np.set_printoptions(suppress=True) # 数值不以科学计数法输出

print(result.x)

print(result.fun)

输出为

[-0.00000003 -0.00000005]
10.000000000000009

去除误差，答案是正确的，也就是在 x=0, y=0，函数 f(x,y) 最小为 10.

如果预设 x0 = np.array([0,0]), 结果就是

[0. 0.]
10.0

1 2	[0. 0.] 10.0

初始预计值还是要有所讲究

求正弦函数的最小值

正弦函数 $f(x)=sin(x)+2$ 是一个周期波，它有无数的波峰和波谷。我们的目的就是要找到猜测试值附近的波谷所在。上代码

np.set_printoptions(suppress=True)

def objective(x):
    return np.sin(x) + 2

x0 = np.array([0])
result = opt.minimize(objective, x0)
print(result.x)
print(result.fun)

print('---')

x0 = np.array([5])
result = opt.minimize(objective, x0)
print(result.x)
print(result.fun)

np.set_printoptions(suppress=True)

def objective(x):

return np.sin(x) + 2

x0 = np.array([0])

result = opt.minimize(objective, x0)

print(result.x)

print(result.fun)

print('---')

x0 = np.array([5])

result = opt.minimize(objective, x0)

print(result.x)

print(result.fun)

分别得到两组结果

[-1.57079633]
1.0
---
[4.71238898]
1.0

最小值都是一样的，所以都找到了各自猜测值附近的波谷位置

补上它的平面图

optimize.minimize 还可以实现曲线拟合问题，比如说找到相近的变化趋势，指纹/人脸检测等，当然这方面有更专业的工具箱

关于背包问题

背包问题是一个整数的线性规划，scipy.optimize.linprog 可用来解决线性规划，但不适于整数的线性规划，因为背包中能装的物品数量是自然数。更适于解决整数的线程规划应该是 Pulp(PuLP is an linear and mixed integer programming modeler) 或 Pyomo.

什么是背包问题:

有 n 种物品，物品 j 的重量和价值分别为 $w_{j}$ 和 $v_{j}$ ，如果背包的最大容量限制是 m ，怎么样选择放入背包的物品以使得背包的总价值最大？

它是一个组合优化问题，设 $x_{j}$ 表示装入背包的第 j 个物品的数量, 那么目标函数和约束条件是：

目标函数： $max\sum_{j=1}^{n}v_{j}x_{j}$
约束条件: $\begin{cases}\sum_{j=1}^n w_j x_j \leq b \\x_j \in \mathbb{N}\end{cases}$

下面是具体问题的求解

from pulp import LpMaximize, LpProblem, LpVariable, lpSum, PULP_CBC_CMD

weight_values = [(2, 3), (3, 4), (4, 5), (5, 6)]  # 每个物品重量及价值
capacity = 5  # 背包的最大承重

problem = LpProblem(name="knapsack-problem", sense=LpMaximize) # 目标函数最大化

# 决策, 1 为选择
x = [LpVariable(f"x{i}", cat="Binary") for i in range(len(weight_values))]

# 目标与约束
problem.objective = lpSum(v[1] * x[i] for i, v in enumerate(weight_values))
problem.addConstraint(lpSum(v[0] * x[i] for i, v in enumerate(weight_values)) <= capacity)

status = problem.solve(PULP_CBC_CMD(msg=False))

print("Selected items:")
for i, v in enumerate(weight_values):
    if x[i].value() == 1:
        print(f" - Item {i + 1} (weight={v[0]}, value={v[1]})")

print("Total value: %d" % problem.objective.value())

from pulp import LpMaximize, LpProblem, LpVariable, lpSum, PULP_CBC_CMD

weight_values = [(2, 3), (3, 4), (4, 5), (5, 6)] # 每个物品重量及价值

capacity = 5 # 背包的最大承重

problem = LpProblem(name="knapsack-problem", sense=LpMaximize) # 目标函数最大化

# 决策, 1 为选择

x = [LpVariable(f"x{i}", cat="Binary") for i in range(len(weight_values))]

# 目标与约束

problem.objective = lpSum(v[1] * x[i] for i, v in enumerate(weight_values))

problem.addConstraint(lpSum(v[0] * x[i] for i, v in enumerate(weight_values)) <= capacity)

status = problem.solve(PULP_CBC_CMD(msg=False))

print("Selected items:")

for i, v in enumerate(weight_values):

if x[i].value() == 1:

print(f" - Item {i + 1} (weight={v[0]}, value={v[1]})")

print("Total value: %d" % problem.objective.value())

解出问题后，结果为

Selected items:
- Item 1 (weight=2, value=3)
- Item 2 (weight=3, value=4)
Total value: 7

最后包的总重为: 5

这个问题还可以更复杂一些，比如说同一件物品可以拿多个进行组合。

链接：1. 最优化学习笔记

本文链接 https://yanbin.blog/scipy-optimize-minimize/, 来自隔叶黄莺 Yanbin Blog

seetimee on 体验 Python FastAPI 的并发能力及线, 进程模型感谢
Yanbin on Mockito 3.4.0 开始可 Mock 静态方法有一个补救，新写了一篇 https://yanbin.blog/mockito-mock-static-method-in-multiple...
Yanbin on 升级到 Spring Boot 3 后 javax.inject.Named 不可用怎么，被抄袭了！算是被机器翻译引用的？
oleksi on 升级到 Spring Boot 3 后 javax.inject.Named 不可用https://www.springcloud.io/post/2023-02/springboot-3-javax-inject-name...
Zhe on 体验 Python FastAPI 的并发能力及线, 进程模型你在trio下得用 trio.sleep()，调用asyncio.sleep() 就会报错。