是否有numpy/scipy函数来计算出站罚款？

Question

在我的极小化问题中，所有有界极小化方法，如“L-BFGS-B”、“TNC”都不收敛，但“Nelder-Mead”收敛得非常好。因此，我更喜欢使用“Nelder Mead”，并修改最小化函数，如下所示：

def outbound_penalty(x, bounds):
    o1 = (bounds[:, 0]-x).max()
    o2 = (x-bounds[:, 1]).max()
    outbound = max(o1, o2, 0)
    rez = 100500*outbound

def bounded_fun(x, bounds):
    return fun(x) + outbound_penalty(x, bounds)

x是numpy数组形状（4），边界有形状（2，4），边界[0]是底部边界，边界[1]-顶部边界。 出乎意料的是，它没有我预料的那么快。在4*10^6次呼叫时，CPU占用自己的时间为40秒。 我当然记得。但我不得不问。numpy/scipy中是否有一些非常优化的函数，可以用来构造出站惩罚

sss = np.zeros((2, 1000))
sss[0] = np.random.uniform(-100, 300, 1000)
sss[1] = np.random.uniform(-100, 300, 1000)
smpls = sss.T

bnd = np.array([[0, 100+np.random.randint(100)], [0, 100+np.random.randint(100)]])
np_bounds = np.array(bnd)

def outbound_penalty(x, bs):
    o1 = (bs[:, 0] - x).max()
    o2 = (x - bs[:, 1]).max()
    outbound = max(o1, o2, 0)
    return 1000000 * outbound

def outbound_penalty_fast(x, bs):
    o1 = (bs[:, 0, None] - x).max(axis=0)
    o2 = (x - bs[:, 1, None]).max(axis=0)
    outbound = np.clip(np.maximum(o1, o2), a_max=None, a_min=0)
    return 1000000 * outbound

%timeit [outbound_penalty(x, np_bounds) for x in smpls]
22.6 ms ± 198 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit [outbound_penalty_fast(x, np_bounds) for x in smpls]
68.6 ms ± 1.62 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit outbound_penalty(smpls[0], np_bounds)
22.5 µs ± 109 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit outbound_penalty_fast(smpls[0], np_bounds)
68.2 µs ± 1.39 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

此版本的出站\u惩罚\u fast返回单浮点数，正如调用者所期望的：

def outbound_penalty_fast(x, bs):
    o1 = (bs[:, 0, None] - x).max(axis=0)
    o2 = (x - bs[:, 1, None]).max(axis=0)
    outbound = np.clip(np.maximum(o1, o2), a_max=None, a_min=0)
    rez = 1000000 * outbound.max()
    return rez

Answer 1

您可以使用广播一次性执行outbound elementwise函数调用，当然也可以使用np.max()，而不是在for循环中比较y和mx：

import numpy as np


def outbound_penalty(x, bs):
    o1 = (bs[:, 0] - x).max()
    o2 = (x - bs[:, 1]).max()
    outbound = max(o1, o2, 0)
    return 1000000 * outbound


def outbound_penalty_fast(x, bs):
    o1 = (bs[:, 0, None] - x).max(axis=0)
    o2 = (x - bs[:, 1, None]).max(axis=0)
    outbound = np.clip(np.maximum(o1, o2), a_max=None, a_min=0)
    return 1000000 * outbound


bnd = np.random.randint(100, 200, size=(2, 2))
bnd[:, 0] = 0
sss = np.random.uniform(-100, 300, size=(2, 1000))

%timeit np.max(np.array([outbound_penalty(x, bnd) for x in sss.T]))
# 9.44 ms ± 166 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit np.max(outbound_penalty_fast(sss, bnd))
# 38.1 µs ± 1.33 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

mx = np.max(np.array([outbound_penalty(x, bnd) for x in sss.T]))
mx_fast = np.max(outbound_penalty_fast(sss, bnd))

mx == mx_fast  # True