是否有一个 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)，bounds 具有形状 (2, 4)，bounds[0] 是底部边框，bounds[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)

此版本的 outbound_penalty_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

您可以使用广播以一次性执行出站元素函数调用，当然当然使用np.max()而不是将y到mx for loop：

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