Python有效地找到多个多项式的局部最大值/最小值

我正在寻找一种有效的方法来查找多个（> 100万）但独立的四阶多项式的局部最小值给定/指定范围/边界.

我有两个要求：

R1：即使对于 100 万个不同的多项式方程也有效

R2：局部最小值精确到 0.01（即 2dp）

这是我使用创建的一些代码scipy。没关系，但我想知道是否有其他更好的软件包可以执行此类任务在我进行并行编程之前。

为了说明我的问题，让我们首先从一个多项式开始：

下面我试图在 (-5, 5) 范围内找到 4x^4 + 6x^3 + 3x^2 + x + 5 的局部最小值。

在我的笔记本电脑上，大约需要 2 毫秒才能找到本地最小值（约为 -0.72770502）。

对于一个多项式来说时间还可以，但我想要更快的速度，因为我需要定期执行此操作超过 100 万次。

from scipy import optimize
import numpy as np

# Define a objective and gradient function for 4th order polynomial
# x is the value to be evaluated
# par is a numpy array of len 5 that specifies the polynomial coefficients.
def obj_grad_fun_custom(x,par):
    obj = (np.array([x**4,x**3,x**2,x**1,1]) * par).sum()
    grad = (np.array([4*x**3,3*x**2,2*x,1]) * par[:-1]).sum()
    return obj, grad

# Try minimise an example polynomial of 4x^4 + 6x^3 + 3x^2 + x + 5
# with contrainted bound
res = optimize.minimize(
    fun = obj_grad_fun_custom,
    x0 = 0,
    args=(np.array([4,6,3,1,5])), # polynomial coefficients
    jac=True ,
    bounds=[(-2, 10)],
    tol=1e-10)
print(res.x)

# Timing (this takes about 2 ms for me)
%timeit optimize.minimize(fun = obj_grad_fun_custom, x0 = 0, args=(np.array([4,6,3,1,5])), jac=True, bounds=[(-5, 5)], tol=1e-10)

下面是我计划对 100 万个不同的 4 阶多项式进行常规操作，我希望在本地最小化这些多项式。希望有人能给我指出一个比scipy。或者还有其他替代方法吗？谢谢！

# Multiple polynomials
result = [] # saving the local minima
poly_sim_no = 1000000 #ideally 1 million or even more
np.random.seed(0)
par_set = np.random.choice(np.arange(10), size=(poly_sim_no, 5), replace=True) #generate some order 4 polynomial coefficients 

for a in par_set:
    res = optimize.minimize(obj_grad_fun_custom, 0,args=(a), jac=True ,bounds=[(-5, 5)], tol=1e-10)
    result.append(res.x)

print(result)

由于您要查找多项式的最小值，因此您可以利用以下事实：获取多项式的导数很容易，并且有许多很好的算法可以查找多项式的根。

它的工作原理如下：

1. 首先，求导数。所有最小值点的导数都为零。
2. 寻找那些零，也就是找到导数的根。
3. 一旦我们有了候选者名单，请检查解决方案是否真实。
4. 检查解决方案是否在您设置的范围内。 （我不知道你添加边界是因为你实际上想要边界，还是为了让它运行得更快。如果是后者，请随意删除这一步。）
5. 实际上用多项式评估候选者并找到最小的一项。

这是代码：

import numpy as np
from numpy.polynomial import Polynomial

def find_extrema(poly, bounds):
    deriv = poly.deriv()
    extrema = deriv.roots()
    # Filter out complex roots
    extrema = extrema[np.isreal(extrema)]
    # Get real part of root
    extrema = np.real(extrema)
    # Apply bounds check
    lb, ub = bounds
    extrema = extrema[(lb <= extrema) & (extrema <= ub)]
    return extrema

def find_minimum(poly, bounds):
    extrema = find_extrema(poly, bounds)
    # Note: initially I tried taking the 2nd derivative to filter out local maxima.
    # This ended up being slower than just evaluating the function.

    # Either bound could end up being the minimum. Check those too.
    extrema = np.concatenate((extrema, bounds))
    # Check every candidate by evaluating the polynomial at each possible minimum,
    # and picking the minimum.
    value_at_extrema = poly(extrema)
    minimum_index = np.argmin(value_at_extrema)
    return extrema[minimum_index]

# Warning: polynomial expects coeffients in the opposite order that you use.
poly = Polynomial([5,1,3,6,4]) 
print(find_minimum(poly, (-5, 5)))

在我的计算机上这需要 162 微秒，比 scipy.optimize 解决方案快大约 6 倍。 （问题中显示的解决方案在我的计算机上需要 1.12 毫秒。）

编辑：更快的替代方案

这是一种更快的方法。然而，它放弃了边界检查，使用了已弃用的 API，并且通常更难以阅读。

p = np.poly1d([4,6,3,1,5])  # Note: polynomials are opposite order of before

def find_minimum2(poly):
    roots = np.real(np.roots(poly.deriv()))
    return roots[np.argmin(poly(roots))]

print(find_minimum2(p))

其计时速度为 110 微秒，比原始速度快大约 10 倍。