可变数量因子的Numpy矢量化求和

Question

我目前正在计算一个包含索引求和的函数。索引在0到T的整数部分之间;理想情况下，我希望能够快速计算T的几个值的总和。 在现实生活中，T的大部分值都很小，但是一小部分可能比T大一个或两个数量级。平均。

我现在在做的是：
1）我定义了向量T，例如， （我的现实生活中的数据具有条目的数目大得多，这是只给一个想法）：

import numpy as np 
T = np.random.exponential(5, 10) 

2）I创建包含0和INT（T之间的因子的矩阵），然后零：

n = int(T.max()) 
j = ((np.arange(n) < T[:,np.newaxis])*np.arange(1,n+1)).astype(int).transpose() 
print(j) 

[[1 1 1 1 1 1 1 1 1 1]
[2 0 2 2 2 0 2 0 2 2]
[0 0 3 0 3 0 3 0 3 3 ]
[0 0 4 0 4 0 0 0 4 4]
[0 0 5 0 5 0 0 0 5 5]
[0 0 6 0 6 0 0 0 6 6]
[0 0 7 0 7 0 0 0 0 7]
[0 0 8 0 8 0 0 0 0 8]
[0 0 9 0 9 0 0 0 0 9]
[0 0 0 0 10 0 0 0 0 10]
[0 0 0 0 11 0 0 0 0 0]
[0 0 0 0 12 0 0 0 0 0]]

3）I生成的总和的单个元件，使用掩模，以避免将所述函数应用于为0的元素：

A = np.log(1 + (1 + j) * 5)* (j>0) 

4）我沿着列总结：

A.sum(axis=0) 

获得： 阵列（[5.170484，2.39789527，29.96464821，5.170484， 42.29052851，2.39789527，8.21500643，2.39789527， 18.49060911，33.9899999]）

是否有最快/更好的矢量化方法？由于大量的零不会导致总和，所以我觉得它很慢，但是因为我是NumPy的初学者，所以我找不出一个更好的写作方法。

编辑：在我的实际问题中，应用于j的函数还取决于第二个参数tau（在相同大小的T的向量中）。所以每列中包含的项目都不相同。

Answer 1

看着你j，为每列有大量的学生都从1到N，其中N正在根据各T元素决定。然后，沿着每一列求和，这与总计相同，直到N，因为无论如何其余的元素都是零。这些总和值可以用np.cumsum和N这些值来计算，基本上j中每列的限制值可以直接从T计算出来。然后将这些N值用作索引来索引cumsum-ed值，以便为我们提供最终输出。

由于cumsum是唯一的计算，所以在一维数组上，与沿着每列的二维数组上的原始方法进行的总和相比，这应该是相当快的和高效的内存。因此，我们有像这样一个量化的方法 -

n = int(T.max()) 
vals = (np.log(1 + (1 + np.arange(1,n+1)) * 5)).cumsum() 
out = vals[(T.astype(int)).clip(max=n-1)] 

在内存使用方面，我们产生三个变量 -

n : Scalar 
vals : 1D array of n elements 
out : 1D array of T.size elements (this is the output anyway) 

运行时测试和验证输出 -

In [5]: def original_app(T): 
    ...:  n = int(T.max()) 
    ...:  j = ((np.arange(n) < T[:,None])*np.arange(1,n+1)).astype(int).transpose() 
    ...:  A = np.log(1 + (1 + j) * 5)* (j>0) 
    ...:  return A.sum(axis=0) 
    ...: 
    ...: def vectorized_app(T): 
    ...:  n = int(T.max()) 
    ...:  vals = (np.log(1 + (1 + np.arange(1,n+1)) * 5)).cumsum() 
    ...:  return vals[(T.astype(int)).clip(max=n-1)] 
    ...: 

In [6]: # Input array 
    ...: T = np.random.exponential(5, 10000) 

In [7]: %timeit original_app(T) 
100 loops, best of 3: 9.62 ms per loop 

In [8]: %timeit vectorized_app(T) 
10000 loops, best of 3: 50.1 µs per loop 

In [9]: np.allclose(original_app(T),vectorized_app(T)) # Verify outputs 
Out[9]: True