问题描述
我有以下多项式方程,我想为其找到局部最小值和最大值。
我将函数定义如下。它使用 flatten 函数来展平嵌套列表,我会将它包含在内以进行测试(在此处找到它 http://rightfootin.blogspot.com/2006/09/more-on-python-flatten.html)
扁平化列表
from itertools import combinations
import math
def flatten(l,ltypes=(list,tuple)):
ltype = type(l)
l = list(l)
i = 0
while i < len(l):
while isinstance(l[i],ltypes):
if not l[i]:
l.pop(i)
i -= 1
break
else:
l[i:i + 1] = l[i]
i += 1
return ltype(l)
我的多项式
def poly(coefficients,factors):
#quadratic terms
constant = 1
singles = factors
products = [math.prod(c) for c in combinations(factors,2)]
squares = [f**2 for f in factors]
sequence = flatten([constant,singles,products,squares])
z = sum([math.prod(i) for i in zip(coefficients,sequence)])
return z
它接受的参数是一个系数列表,例如:
coefs = [12.19764959,-1.8233151,2.50952816,-1.56344375,1.00003828,-1.72128301,-2.54254877,-1.20377309,5.53510616,2.94755653,4.83759279,-0.85507208,-0.48007208,-3.70507208,-0.27007208]
以及因子或变量值的列表:
factors = [0.4714,0.4714,-0.4714,0.4714]
将这些插入并计算多项式的结果。我这样写的原因是因为变量(因子)的数量从适合到适合变化,所以我想保持它的灵活性。我现在想在函数达到其最大值和最小值的某个范围内(假设在 -1 和 1 之间)找到“因子”值的组合。如果函数是“硬编码”的,我可以使用 scipy.optimize,但我不知道如何让它按原样工作。
另一种选择是蛮力网格搜索(我目前使用),但是一旦您有超过 2 个变量,它就会非常慢,尤其是小步长。可能没有真正的最小值/最大值,其中斜率 == 0 在界限内,但只要我能得到最大值和最小值就可以了。
解决方法
好的,我想通了。这是两件非常愚蠢的事情:
-
函数中参数的顺序必须颠倒,以便第一个参数(我想要优化的参数)是“因子”或 X 值,然后是系数。这样就可以将相同大小的数组用作 X0,并将系数用作参数。
-
这还不够,因为如果输入的是数组,该函数将返回一个数组。我只是在函数本身中添加了一个 factor = list(factors) 以使其成为正确的形状。
新功能:
def poly(factors,coefficients):
factors = list(factors)
#quadratic terms
constant = 1
singles = factors
products = [math.prod(c) for c in combinations(factors,2)]
squares = [f**2 for f in factors]
sequence = flatten([constant,singles,products,squares])
z = sum([math.prod(i) for i in zip(coefficients,sequence)])
return z
还有优化:
coefs = [4.08050532,-0.47042713,-0.08200181,-0.54184481,-0.18515675,-0.96751856,-1.10814625,-1.7831592,5.2763512,2.83505438,4.7082153,0.22988773,1.06488773,-0.70011227,1.42988773]
x0 = [0.1,0.1,0.1]
minimize(poly,x0 = x0,args = coefs,bounds = ((-1,1),(-1,1)))
哪个返回:
fun: -1.6736636102536673
hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
jac: array([-2.10611305e-01,2.19138777e+00,-8.16990766e+00,-1.11022302e-07])
message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
nfev: 85
nit: 12
njev: 17
status: 0
success: True
x: array([1.,-1.,1.,0.03327357])