标签: vectorization

背景

我的问题是由简单的观察推动的,这有点破坏了经验丰富的MATLAB用户经常持有/做出的信念/假设:

• 当涉及内置函数和基本语言特性(如索引向量和矩阵)时,MATLAB得到了很好的优化.
• MATLAB中的循环很慢(尽管有JIT),如果算法可以以原生的"矢量化"方式表示,通常应该避免.

底线:核心MATLAB功能是高效的,并且尝试使用MATLAB代码超越它是很困难的,如果不是不可能的话.

研究向量索引的性能

下面显示的示例代码是基本的:我为所有向量条目指定标量值.首先,我分配一个空向量x:

tic; x = zeros(1e8,1); toc
Elapsed time is 0.260525 seconds.

有x,我想它所有的项目设置为相同的值.在实践中,你会采用不同的方式,例如,x = value*ones(1e8,1)但这里的重点是研究向量索引的性能.最简单的方法是写:

tic; x(:) = 1; toc
Elapsed time is 0.094316 seconds.

我们称之为方法1(从赋值的值x).它看起来非常快(至少比内存分配更快).因为我在这里做的唯一事情就是对内存进行操作,我可以通过计算获得的有效内存带宽并将其与计算机的硬件内存带宽进行比较来估算此代码的效率:

eff_bandwidth = numel(x) * 8 bytes per double * 2 / time

在上面,我乘以2因为除非使用SSE流,否则在存储器中设置值要求向量既可以从存储器读取也可以写入存储器.在上面的例子中:

eff_bandwidth(1) = 1e8*8*2/0.094316 = 17 Gb/s

我的计算机的STREAM基准内存带宽约为17.9 Gb/s,所以确实如此 - 在这种情况下,MATLAB可以提供接近峰值的性能!到现在为止还挺好.

Method 1 is suitable if you want to set all …

为什么,在硬件执行操作的最低级别和所涉及的一般底层操作(即:运行代码时所有编程语言的实际实现的一般性),矢量化通常比循环更快？

当使用矢量化时,计算机在循环时做了什么(我说的是计算机执行的实际计算,而不是程序员编写的计算),或者它有什么不同的做法？

我一直无法说服自己为什么差异应该如此重要.我可能会说服矢量化代码在某处削减一些循环开销,但计算机仍然必须执行相同数量的操作,不是吗？例如,如果我们将大小为N的向量乘以标量,我们将使用N次乘法执行任一方式,不是吗？

i=np.arange(1,4,dtype=np.int)
a=np.arange(9).reshape(3,3)

和

a
>>>array([[0, 1, 2],
          [3, 4, 5],
          [6, 7, 8]])
a[:,0:1]
>>>array([[0],
          [3],
          [6]])
a[:,0:2]
>>>array([[0, 1],
          [3, 4],
          [6, 7]])
a[:,0:3]
>>>array([[0, 1, 2],
          [3, 4, 5],
          [6, 7, 8]])

现在我想对数组进行矢量化以将它们一起打印出来.我试试

a[:,0:i]

要么

a[:,0:i[:,None]]

它给出了TypeError:只有整数标量数组可以转换为标量索引

我想知道list_1中包含哪些元素list_2。我需要输出作为布尔值的有序列表。但我想避免for循环，因为两个列表都有超过 200 万个元素。

这就是我所拥有的并且它有效，但它太慢了：

list_1 = [0,0,1,2,0,0]
list_2 = [1,2,3,4,5,6]

booleans = []
for i in list_1:
   booleans.append(i in list_2)

# booleans = [False, False, True, True, False, False]

我可以拆分列表并使用多线程，但如果可能的话，我更喜欢更简单的解决方案。我知道像 sum() 这样的一些函数使用向量运算。我正在寻找类似的东西。

如何让我的代码更加高效？

假设我有一个NxN矩阵A,一个由数字1:N的子集和一个值K组成的索引向量V,我想这样做:

 for i = V
     A(i,i) = K
 end

有没有办法在一个带矢量化的声明中做到这一点？

例如A(某事)= K.

该声明A(V,V) = K不起作用,它分配非对角元素,这不是我想要的.例如:

>> A = zeros(5);
>> V = [1 3 4];
>> A(V,V) = 1

A =

 1     0     1     1     0
 0     0     0     0     0
 1     0     1     1     0
 1     0     1     1     0
 0     0     0     0     0

x <- seq(0.1,10,0.1)
y <- if (x < 5) 1 else 2

我希望if能够对每一个案例进行操作,而不是对整个向量进行操作.我需要改变什么？

之前在比较bsxfun和repmat表现之间提出的问题很少.

• 其中一个是:Matlab - bsxfun no longer faster than repmat?.这个尝试研究repmat和之间的性能比较bsxfun,特别是从输入数组本身的列中执行减去输入数组的平均值,因此只会探索与其等价物相对的@minus部分.bsxfunrepmat
• 另一个是:In Matlab, when is it optimal to use bsxfun?.那个试图通过沿列的平均值进行相同的减法操作,也没有扩展到其他内置操作.

在这篇文章中,我试图调查两者之间的性能数据bsxfun并repmat覆盖所有bsxfun内置函数,从而为它提供更广泛的视角,因为这两者都提供了良好的矢量化解决方案.

具体来说,我对这篇文章的问题是:

1. 各种内置操作如何bsxfun对repmat等效执行？bsxfun支持浮点运算一样@plus,@minus,@times等,并且还像关系和逻辑运算@ge,@and等等.所以,是否有特定的内置插件,会给我明显的加速与bsxfun比使用他们的repmat等价物？

2. 罗兰在她blog post已经基准repmat对bsxfun具有定时@() A - repmat(mean(A),size(A,1),1)对@() bsxfun(@minus,A,mean(A))分别.如果我需要涵盖所有内置函数的基准测试,我可以使用一些其他可用于浮点,关系和逻辑运算的比较模型吗？

我很感激在寻找和理解pythonic方法的一些帮助,以优化嵌套for循环中的以下数组操作:

def _func(a, b, radius):
    "Return 0 if a>b, otherwise return 1"
    if distance.euclidean(a, b) < radius:
        return 1
    else:
        return 0

def _make_mask(volume, roi, radius):
    mask = numpy.zeros(volume.shape)
    for x in range(volume.shape[0]):
        for y in range(volume.shape[1]):
            for z in range(volume.shape[2]):
                mask[x, y, z] = _func((x, y, z), roi, radius)
    return mask

其中volume.shape(182,218,200)和roi.shape(3,)都是ndarray类型; 而且radius是一个int

我有一个相对稀疏的数组,我想遍历每一行并只调整非零元素.

示例输入:

[2,3,1,0]
[0,0,2,1]

示例输出:

[2,1,3,0]
[0,0,1,2]

注意零如何改变位置.

要对每行中的所有元素(包括零)进行洗牌,我可以这样做:

for i in range(len(X)):
    np.random.shuffle(X[i, :])

我试图做的是这样的:

for i in range(len(X)):
    np.random.shuffle(X[i, np.nonzero(X[i, :])])

但它没有效果.我注意到返回类型与X[i, np.nonzero(X[i, :])]原因不同X[i, :].

In[30]: X[i, np.nonzero(X[i, :])]
Out[30]: array([[23,  5, 29, 11, 17]])

In[31]: X[i, :]
Out[31]: array([23,  5, 29, 11, 17])

想象一下,你有一个很长的序列.查找序列全为零的间隔的最有效方法是什么(或者更确切地说,序列降至接近零的值abs(X)<eps):

为简单起见,我们假设以下顺序:

sig = [1 1 0 0 0 0 1 1 1 1 1 0 1 0 0 0 1 1 1 1 1 1 1 1 0 0 1 1 1 0];

我正在尝试获取以下信息:

startIndex   EndIndex    Duration
3            6           4
12           12          1
14           16          3
25           26          2
30           30          1

然后使用这些信息,我们找到持续时间> =到某个指定值(例如3)的间隔,并返回所有这些间隔中的值的索引组合:

indices = [3 4 5 6 14 15 16];

最后一部分与前一个问题有关:

MATLAB:从开始/结束索引列表创建矢量化数组

这是我到目前为止:

sig = [1 1 0 0 0 0 1 …