numpy 向量化多维函数

Question

我在矢量化多维函数时遇到问题。
考虑以下示例：

def _cost(u):
    return u[0] - u[1]

cost = np.vectorize(_cost)

>>> x = np.random.normal(0, 1,(10, 2))
>>> cost(x)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/lucapuggini/MyApps/scientific_python_3_5/lib/python3.5/site-packages/numpy/lib/function_base.py", line 2218, in __call__
    return self._vectorize_call(func=func, args=vargs)
  File "/Users/lucapuggini/MyApps/scientific_python_3_5/lib/python3.5/site-packages/numpy/lib/function_base.py", line 2281, in _vectorize_call
    ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args)
  File "/Users/lucapuggini/MyApps/scientific_python_3_5/lib/python3.5/site-packages/numpy/lib/function_base.py", line 2243, in _get_ufunc_and_otypes
    outputs = func(*inputs)
TypeError: _cost() missing 1 required positional argument: 'v'

背景信息：我在尝试将以下代码（粒子群优化算法）推广到多元数据时遇到了问题：

import numpy as np
import matplotlib.pyplot as plt


def pso(cost, sim, space_dimension, n_particles, left_lim, right_lim, f1=1, f2=1, verbose=False):

    best_scores = np.array([np.inf]*n_particles)
    best_positions = np.zeros(shape=(n_particles, space_dimension))
    particles = np.random.uniform(left_lim, right_lim, (n_particles, space_dimension))
    velocities = np.zeros(shape=(n_particles, space_dimension))

    for i in range(sim):
        particles = particles + velocities
        print(particles)
        scores = cost(particles).ravel()
        better_positions = np.argwhere(scores < best_scores).ravel()
        best_scores[better_positions] = scores[better_positions]
        best_positions[better_positions, :] = particles[better_positions, :]
        g = best_positions[np.argmin(best_scores), :]

        u1 = np.random.uniform(0, f1, (n_particles, 1))
        u2 = np.random.uniform(0, f2, (n_particles, 1))
        velocities = velocities + u1 * (best_positions - particles) + u2 * (g - particles)

        if verbose and i % 50 == 0:
            print('it=', i, ' score=', cost(g))


            x = np.linspace(-5, 20, 1000)
            y = cost(x)

            plt.plot(x, y)
            plt.plot(particles, cost(particles), 'o')
            plt.vlines(g, y.min()-2, y.max())
            plt.show()


    return g, cost(g)




def test_pso_1_dim():

    def _cost(x):
        if 0 < x < 15:
            return np.sin(x)*x 
        else:
            return 15 + np.min([np.abs(x-0), np.abs(x-15)])

    cost = np.vectorize(_cost)

    sim = 100
    space_dimension = 1
    n_particles = 5
    left_lim, right_lim = 0, 15
    f1, f2  = 1, 1

    x, cost_x = pso(cost, sim, space_dimension, n_particles,
                    left_lim, right_lim, f1, f2, verbose=False)

    x0 = 11.0841839
    assert np.abs(x - x0) < 0.01

    return 

如果在这种情况下矢量化不是一个好主意，请告诉我。

Answer 1

正如注释中提到的vectorize：

提供矢量化函数主要是为了方便，而不是为了性能。该实现本质上是一个 for 循环。

因此，虽然通过numpy类型和函数对代码进行矢量化可能是一个好主意，但您可能不应该使用numpy.vectorize.

对于您给出的示例，您cost可以简单有效地计算为在数组上运行的函数numpy：

def cost(x):
    # Create the empty output 
    output = np.empty(x.shape)
    # Select the first group using a boolean array
    group1 = (0 < x) & (x < 15) 
    output[group1] = np.sin(x[group1])*x[group1]
    # Select second group as inverse (logical not) of group1 
    output[~group1] = 15 + np.min(
        [np.abs(x[~group1]-0), np.abs(x[~group1]-15)],
        axis=0)
    return output