【问题标题】：How to vectorize finding max value in numpy array with if statement?如何使用 if 语句矢量化在 numpy 数组中查找最大值？
【发布时间】：2014-01-08 01:06:45
【问题描述】：

我的设置：Python 2.7.4.1、Numpy MKL 1.7.1、Windows 7 x64、WinPython

上下文：

我尝试实现用于求解 SVM 的序列最小优化算法。我使用最大违例对方法。

问题：

在工作集选择过程中，我想找到满足某些条件的元素的梯度最大值及其索引，y[i]*alpha[i]

#y - array of -1 and 1
y=np.array([-1,1,1,1,-1,1])
#alpha- array of floats in range [0,C]
alpha=np.array([0.4,0.1,1.33,0,0.9,0])
#grad - array of floats
grad=np.array([-1,-1,-0.2,-0.4,0.4,0.2])


GMaxI=float('-inf')       
GMax_idx=-1        
n=alpha.shape[0] #usually n=100000
C=4
B=[0,0,C]
for i in xrange(0,n):
    yi=y[i]  #-1 or 1
    alpha_i=alpha[i]

    if (yi * alpha_i< B[yi+1]): # B[-1+1]=0 B[1+1]=C
        if( -yi*grad[i]>=GMaxI):
            GMaxI= -yi*grad[i]
            GMax_idx = i

这个过程被调用了很多次（~50000），分析器显示这是瓶颈。这段代码可以向量化吗？

编辑 1：添加一些小的示例数据

编辑 2：我检查了 hwlau 、 larsmans 和 E 先生提出的解决方案。只有 E 先生提出的解决方案是正确的。下面是所有三个答案的示例代码：

import numpy as np
y=np.array([   -1,  -1,   -1,  -1,   -1,   -1,   -1,   -1])
alpha=np.array([0,  0.9,  0.4, 0.1, 1.33,    0,  0.9,    0])
grad=np.array([-3, -0.5,   -1,  -1, -0.2,   -4, -0.4, -0.3])
C=4
B=np.array([0,0,C])

#hwlau - wrong index and value
filter = (y*alpha < C*0.5*(y+1)).astype('float')
GMax_idx = (filter*(-y*grad)).argmax()
GMax = -y[GMax_idx]*grad[GMax_idx]

print GMax_idx,GMax

#larsmans - wrong index
neg_y_grad = (-y * grad)[y * alpha < B[y + 1]]
GMaxI = np.max(neg_y_grad)
GMax_ind = np.argmax(neg_y_grad)

print GMax_ind,GMaxI

#Mr E - correct result
BY = np.take(B, y+1)
valid_mask = (y * alpha < BY)
values = -y * grad
values[~valid_mask] = np.min(values) - 1.0

GMaxI = values.max()
GMax_idx = values.argmax()
print GMax_idx,GMaxI

Output (GMax_idx, GMaxI)
0 -3.0 
3 -0.2
4 -0.2

结论

检查所有解决方案后，最快的（2x-6x）是@ali_m 提出的解决方案。但是它需要安装一些 python 包：numba 及其所有先决条件。

我在使用 numba 和类方法时遇到了一些麻烦，所以我创建了使用 numba 自动处理的全局函数，我的解决方案如下所示：

from numba import autojit  

@autojit
def FindMaxMinGrad(A,B,alpha,grad,y):
    '''
    Finds i,j indices with maximal violatin pair scheme
    A,B - 3 dim arrays, contains bounds A=[-C,0,0], B=[0,0,C]
    alpha - array like, contains alpha coeficients
    grad - array like, gradient
    y - array like, labels
    '''
    GMaxI=-100000
    GMaxJ=-100000

    GMax_idx=-1
    GMin_idx=-1

    for i in range(0,alpha.shape[0]):

        if (y[i] * alpha[i]< B[y[i]+1]):
            if( -y[i]*grad[i]>GMaxI):
                GMaxI= -y[i]*grad[i]
                GMax_idx = i

        if (y[i] * alpha[i]> A[y[i]+1]):
            if( y[i]*grad[i]>GMaxJ):
                GMaxJ= y[i]*grad[i]
                GMin_idx = i

    return (GMaxI,GMaxJ,GMax_idx,GMin_idx)  

class SVM(object):

    def working_set(self,....):
        FindMaxMinGrad(.....)

【问题讨论】：

您可以使用列表理解开始 - 但我怀疑它会使其执行得更快（如果有更快的话）。也许尝试查看内置的 Numpy 函数来完成一些工作。
你能添加一些正确形状和dtype的输入数据作为Python代码吗？例如您可以手动输入n 的小值，也可以使用 NumPy 的随机模块生成正确形状的随机数据。使我们可以直接运行此代码而无需猜测输入和输出值。那么你很快就会得到一个有用的答案。
我已经编辑了我的问题并添加了一些数据。
确实，我在argmax 中有一个错误。解决了这个问题；仍然只有三行代码。

标签： python arrays numpy max vectorization

【解决方案1】：

如果您使用 numba 对使用嵌套循环的原始代码进行 JIT 编译，您可能会比普通矢量化做得更好。

import numpy as np
from numba import autojit

@autojit
def jit_max_grad(y, alpha, grad, B):
    maxgrad = -inf
    maxind = -1
    for ii in xrange(alpha.shape[0]):
        if (y[ii] * alpha[ii] < B[y[ii] + 1]):
            g = -y[ii] * grad[ii]
            if g >= maxgrad:
                maxgrad = g
                maxind = ii
    return maxind, maxgrad

为了比较，这里是 E 先生的矢量化版本：

def mr_e_max_grad(y, alpha, grad, B):

    BY = np.take(B, y+1)
    valid_mask = (y * alpha < BY)
    values = -y * grad
    values[~valid_mask] = np.min(values) - 1.0
    GMaxI = values.max()
    GMax_idx = values.argmax()
    return GMax_idx, GMaxI

时间：

y = np.array([   -1,  -1,   -1,  -1,   -1,   -1,   -1,   -1])
alpha = np.array([0,  0.9,  0.4, 0.1, 1.33,    0,  0.9,    0])
grad = np.array([-3, -0.5,   -1,  -1, -0.2,   -4, -0.4, -0.3])
C = 4
B = np.array([0,0,C])

%timeit mr_e_max_grad(y, alpha, grad, B)
# 100000 loops, best of 3: 19.1 µs per loop

%timeit jit_max_grad(y, alpha, grad, B)
# 1000000 loops, best of 3: 1.07 µs per loop

更新：如果您想查看更大数组上的时序，可以轻松定义一个函数，根据您在问题中的描述生成半真实的假数据：

def make_fake(n, C=4):
    y = np.random.choice((-1, 1), n)
    alpha = np.random.rand(n) * C
    grad = np.random.randn(n)
    B = np.array([0,0,C])
    return y, alpha, grad, B

%%timeit y, alpha, grad, B = make_fake(100000, 4)
mr_e_max_grad(y, alpha, grad, B)
# 1000 loops, best of 3: 1.83 ms per loop

%%timeit y, alpha, grad, B = make_fake(100000, 4)
jit_max_grad(y, alpha, grad, B)
# 1000 loops, best of 3: 471 µs per loop

【讨论】：

谢谢，非常有用的建议。我想知道当数组变大时时间会是什么样子。
实际上这个函数是一个类'SVM'中的一个方法。我尝试使用您的方法，但装饰器 @autojit 不起作用。我对 numba 不是很熟悉，也不知道如何为此类功能添加特定的装饰器。你展示的时间非常令人鼓舞。另外在同一个循环中我必须找到最小值，所以我认为这个解决方案是最好的。
看看这个：numba.pydata.org/numba-doc/0.6/doc/extension_types.html

【解决方案2】：

我认为这是一个完全矢量化的版本

import numpy as np

#y - array of -1 and 1
y=np.array([-1,1,1,1,-1,1])
#alpha- array of floats in range [0,C]
alpha=np.array([0.4,0.1,1.33,0,0.9,0])
#grad - array of floats
grad=np.array([-1,-1,-0.2,-0.4,0.4,0.2])


BY = np.take(B, y+1)
valid_mask = (y * alpha < BY)
values = -yi * grad
values[~valid_mask] = np.min(values) - 1.0

GMaxI = values.max()
GMax_idx = values.argmax()

【讨论】：

感谢您的时间和精力，我会在我的上下文中检查并让您知道计时结果。

【解决方案3】：

给你：

y=np.array([-1,1,1,1,-1,1])
alpha=np.array([0.4,0.1,1.33,0,0.9,0])
grad=np.array([-1,-1,-0.2,-0.4,0.4,0.2])
C=4

filter = (y*alpha < C*0.5*(y+1)).astype('float')
GMax_idx = (filter*(-y*grad)).argmax()
GMax = -y[GMax_idx]*grad[GMax_idx]

没有尝试过基准测试，但它是纯数字和矢量化的，所以应该很快。

【讨论】：

不错且紧凑的解决方案，我将检查所有三个答案，并在我的上下文中让您了解计时结果。

【解决方案4】：

如果将B 从list 更改为NumPy 数组，您至少可以将yi * alpha_i< B[yi+1] 向量化并将循环向内推：

GMaxI = float('-inf')      
GMax_idx = -1      
for i in np.where(y * alpha < B[y + 1])[0]:
    if -y[i] * grad[i] >= GMaxI:
        GMaxI= -y[i] * grad[i]
        GMax_idx = i

这应该会节省一些时间。接下来，可以矢量化-y[i] * grad[i]：

GMaxI = float('-inf')
GMax_idx = -1
neg_y_grad = -y * grad
for i in np.where(y * alpha < B[y + 1])[0]:
    if neg_y_grad[i] >= GMaxI:
        GMaxI= -y[i] * grad[i]
        GMax_idx = i

最后，我们可以通过在-y * grad 上使用max 和argmax 来矢量化整个循环，并由y * alpha < B[y + 1] 过滤：

neg_y_grad = (-y * grad)
GMaxI = np.max(neg_y_grad[y * alpha < B[y + 1]])
GMax_idx = np.where(neg_y_grad == GMaxI)[0][0]

【讨论】：