如何在有界的python优化中找到全局最小值?

Posted

技术标签:

【中文标题】如何在有界的python优化中找到全局最小值?【英文标题】:how to find global minimum in python optimization with bounds? 【发布时间】:2014-03-07 08:43:36 【问题描述】:

我有一个包含 64 个变量的 Python 函数,我尝试在最小化函数中使用 L-BFGS-B 方法对其进行优化,但是该方法对初始猜测有很强的依赖性,未能找到全局最小值.

但我喜欢它为变量设置界限的能力。有没有一种方法/函数可以在变量有边界的同时找到全局最小值?

【问题讨论】:

我怀疑math.stackexchange.com 是更适合提出此类问题的地方。 您能否描述一下您的功能——平滑/渐变/Hessian?如果您可以将其表述为平方和,请参阅scipy-optimize-leastsq-with-bound-constraints。另见scicomp.stackexchange.com/search?q=bfgs。 我在 3D 空间中设计了 8 条贝塞尔曲线,每条曲线有 6 个控制点,要最小化的函数是这些曲线的评价函数,它是 4 个不同参数(长度、半径曲率、接近度、高度顺序)从曲线推导出来。到目前为止,我尝试了 scipy.minimize()、basinhopping,但我仍然无法找到全局最小值 【参考方案1】:

这可以通过scipy.optimize.basinhopping 完成。 Basinhopping 是一种旨在找到目标函数的全局最小值的函数。它使用函数scipy.optimize.minimize 进行重复最小化,并在每次最小化后在坐标空间中采取随机步骤。 Basinhopping 仍然可以通过使用实现边界的最小化器之一(例如 L-BFGS-B)来尊重边界。这里有一些代码展示了如何做到这一点

# an example function with multiple minima
def f(x): return x.dot(x) + sin(np.linalg.norm(x) * np.pi)

# the starting point
x0 = [10., 10.]

# the bounds
xmin = [1., 1.]
xmax = [11., 11.]

# rewrite the bounds in the way required by L-BFGS-B
bounds = [(low, high) for low, high in zip(xmin, xmax)]

# use method L-BFGS-B because the problem is smooth and bounded
minimizer_kwargs = dict(method="L-BFGS-B", bounds=bounds)
res = basinhopping(f, x0, minimizer_kwargs=minimizer_kwargs)
print res

上面的代码适用于一个简单的情况,但如果盆地跳跃随机位移例程将您带到那里,您仍然可能最终进入禁区。幸运的是,可以通过使用关键字 take_step 传递自定义步骤例程来覆盖它

class RandomDisplacementBounds(object):
    """random displacement with bounds"""
    def __init__(self, xmin, xmax, stepsize=0.5):
        self.xmin = xmin
        self.xmax = xmax
        self.stepsize = stepsize

    def __call__(self, x):
        """take a random step but ensure the new position is within the bounds"""
        while True:
            # this could be done in a much more clever way, but it will work for example purposes
            xnew = x + np.random.uniform(-self.stepsize, self.stepsize, np.shape(x))
            if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
                break
        return xnew

# define the new step taking routine and pass it to basinhopping
take_step = RandomDisplacementBounds(xmin, xmax)
result = basinhopping(f, x0, niter=100, minimizer_kwargs=minimizer_kwargs,
                      take_step=take_step)
print result

【讨论】:

无循环的变体:return np.clip( x + random(...), xmin, xmax )。 (在所有维度上使用相同的 step_size,预先缩放,因此边界框大致为正方形。) 你好。我认为您在 call 中需要 self.xmin 和 self.xmax 而不是 xmin 和 xmax?【参考方案2】:

一些关于调试和可视化优化器的常识性建议 关于你的功能:

您的目标函数和约束是否合理? 如果目标函数是一个总和说f() + g(), 为"fx-opt.nptxt"(下)中的所有x单独打印它们; 如果f() 是总和的 99 % 而 g() 是 1 %,则进行调查。

约束:xfinal 中有多少组件 x_i 卡在边界上, x_i &lt;= lo_i&gt;= hi_i ?


您的职能在全球范围内有多坎坷? 使用几个随机起点运行,并保存结果以进行分析/绘图:
title = "%s  n %d  ntermhess %d  nsample %d  seed %d" % (  # all params!
    __file__, n, ntermhess, nsample, seed )
print title
...
np.random.seed(seed)  # for reproducible runs
np.set_printoptions( threshold=100, edgeitems=10, linewidth=100,
        formatter = dict( float = lambda x: "%.3g" % x ))  # float arrays %.3g

lo, hi = bounds.T  # vecs of numbers or +- np.inf
print "lo:", lo
print "hi:", hi

fx = []  # accumulate all the final f, x
for jsample in range(nsample):
        # x0 uniformly random in box lo .. hi --
    x0 = lo + np.random.uniform( size=n ) * (hi - lo)

    x, f, d = fmin_l_bfgs_b( func, x0, approx_grad=1,
                m=ntermhess, factr=factr, pgtol=pgtol )
    print "f: %g  x: %s  x0: %s" % (f, x, x0)
    fx.append( np.r_[ f, x ])

fx = np.array(fx)  # nsample rows, 1 + dim cols
np.savetxt( "fx-opt.nptxt", fx, fmt="%8.3g", header=title )  # to analyze / plot

ffinal = fx[:,0]
xfinal = fx[:,1:]
print "final f values, sorted:", np.sort(ffinal)
jbest = ffinal.argmin()
print "best x:", xfinal[jbest]

如果某些ffinal 值看起来相当不错, 在这些附近尝试更多随机起点- 这肯定比纯随机要好。

如果 x 是曲线或任何真实的,请绘制最好的几个 x0xfinal。 (经验法则是 d 维度中的 nsample ~ 5*d 或 10*d。 太慢了,太多了?减少maxiter / maxeval,减少ftol -- 你不需要ftol 1e-6 来进行这样的探索。)

如果您想要可重现的结果, 那么您必须在title 中列出所有相关参数 并在派生文件和绘图中。 否则,你会问“this 是从哪里来的??”


你的函数在 epsilon 尺度上有多颠簸 ~ 10^-6 ? 近似梯度的方法有时会返回最后的估计, 但如果不是:
from scipy.optimize._numdiff import approx_derivative  # 3-point, much better than
## from scipy.optimize import approx_fprime
for eps in [1e-3, 1e-6]:
    grad = approx_fprime( x, func, epsilon=eps )
    print "approx_fprime eps %g: %s" % (eps, grad)

如果在优化器退出之前梯度估计很差/颠簸, 你不会看到的。 然后你必须保存所有中间的[f, x, approx_fprime] 也看他们;在 python 中很容易——问是否不清楚。

在某些问题区域中,通常会从声称的 xmin 进行备份和重新启动。 例如,如果你在乡间小路上迷路了, 先找一条主路,然后从那里重新开始。


摘要: 不要指望任何黑盒优化器在一个函数上工作 大规模颠簸,或ε级颠簸,或两者兼而有之。 投资于测试脚手架,以及了解优化器在做什么的方式。

【讨论】:

【参考方案3】:

自从您提出这个问题以来,全局优化方面已经取得了一些不错的进展,可能会对您有所帮助。特别是,我会提请您注意 SHGO 算法 (package),它现在也在 scipy.optimize 中作为标准选项之一。但是,如果您真的无法减少它,它可能会与您的搜索空间的维数产生冲突。

你可以尝试一些经典的方法,比如模式搜索,或者PySOT 中的代理方法,benchmark very well 也一样。如果你真的被卡住了,可以考虑 optuna 之类的东西,或者,如果你很绝望,可以考虑 hyperopt。

我目前的选择是 DLIB 和 Nevergrad。两者都相当快。对于没有任何真正依赖关系的选项,也许可以看看 freelunch。

【讨论】:

【参考方案4】:

非常感谢您的详细回复,但由于我对python还很陌生,我不太清楚如何将代码实现到我的程序中,但这是我的优化尝试:

x0=np.array((10, 13, f*2.5, 0.08,    10, f*1.5,  0.06, 20, 
             10, 14, f*2.5, 0.08,    10, f*1.75, 0.07, 20,
             10, 15, f*2.5, 0.08,    10, f*2,    0.08, 20,
             10, 16, f*2.5, 0.08,    10, f*2.25, 0.09, 20,
             10, 17, f*2.5, -0.08,    10, f*2.5, -0.06, 20,
             10, 18, f*2.5, -0.08,    10, f*2.75,-0.07, 20,
             10, 19, f*2.5, -0.08,    10, f*3,   -0.08, 20,
             10, 20, f*2.5, -0.08,    10, f*3.25,-0.09, 20))

# boundary for each variable, each element in this restricts the corresponding element     above
bnds=((1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), )

from scipy.optimize import basinhopping
from scipy.optimize import minimize

merit=a*meritoflength + b*meritofROC + c*meritofproximity +d*(distancetoceiling+distancetofloor)+e*heightorder
minimizer_kwargs = "method": "L-BFGS-B", "bounds": bnds, "tol":1e0
ret = basinhopping(merit_function, x0, minimizer_kwargs=minimizer_kwargs, niter=10, T=0.01)

zoom = ret['x']
res = minimize(merit_function, zoom, method = 'L-BFGS-B', bounds=bnds, tol=1e-5)
print res

评价函数将 x0 与其他一些值结合起来,为 8 条曲线形成 6 个控制点,然后计算它们的长度、曲率半径等。它以这些参数与一些权重的线性组合的形式返回最终的评价。

我使用basinhopping 以低精度找到一些最小值,然后使用minimize 提高最低最小值的精度。

附言我运行的平台是 Enthoght canopy 1.3.0, numpy 1.8.0 scipy 0.13.2 mac 10.8.3

【讨论】:

您的函数通过盆地“跳跃”了多少?对于可能尝试此操作的其他人会很有用。还有一个普遍的问题:控制点是限制在一个大盒子上,还是每个都在自己的 1/6 中?如果是一个大盒子,那么解空间是6! = 比需要大 720 倍。 我如何检查函数跳了多少?至于 6 个控制点,比如 P1、P2、P3、P4、P5、P6。在这些点中,P1、P6 是固定的,P2y、P2z、P5y、p5z 是固定的,因此变量是 P2x、P3x、P3y、P3z、P4x、P4y、P4z、P5x。现在我相信 P3 和 P4 共享同一个大盒子,而 P2、P5 在单个 x 轴上移动 我会提出一个新问题,例如“我如何检查函数的跳跃方式,对于 scipy.optimize.basinhopping”。但是你得到了合理的结果,对吗?

以上是关于如何在有界的python优化中找到全局最小值?的主要内容,如果未能解决你的问题,请参考以下文章

Python 的 hash() 函数返回的最大值/最小值

如何在有约束的 scipy 中使用最小化函数

最优化问题求解方法

Panda DF:在有条件的多列中查找最小值

优化平面任意形状的布局

寻找有界子图之间的最小割集