标签: vectorization

我正在使用 Matlab。

我有一个二元方阵。对于每一行，有一个或多个 1 的条目。我想遍历此矩阵的每一行并返回这些 1 的索引并将它们存储在单元格的条目中。

我想知道是否有一种方法可以在不循环此矩阵的所有行的情况下执行此操作，因为在 Matlab 中 for 循环非常慢。

例如，我的矩阵

M = 0 1 0
    1 0 1
    1 1 1 

然后最终，我想要类似的东西

A = [2]
    [1,3]
    [1,2,3]

A细胞也是如此。

有没有办法在不使用 for 循环的情况下实现这个目标，目的是更快地计算结果？

最初这是mathematica.SE 中出现的一个问题，但由于讨论涉及多种编程语言，我认为最好稍微改一下措辞并将其发布在这里。

简而言之，michalkvasnicka在以下 MATLAB 示例中发现

s = 15000;
tic
% for-loop version
H = zeros(s,s);
for c = 1:s
    for r = 1:s
        H(r,c) = 1/(r+c-1);
    end
end
toc
%Elapsed time is 1.359625 seconds.... For-loop 

tic;
% vectorized version
c = 1:s;
r = c';
HH=1./(r+c-1);
toc
%Elapsed time is 0.047916 seconds.... Vectorized

isequal(H,HH)

矢量化代码段比纯 for 循环代码段快 25 倍以上。虽然我无法访问 MATLAB，因此无法亲自测试示例，但时间1.359625似乎表明它是在普通 PC 上进行测试的，就像我的一样。

但我无法用 Fortran 或 Julia 等其他语言重现时间！（我们知道，他们俩都以数值计算的性能而闻名。好吧，我承认我绝不是 fortran 或 julia 的专家。）

以下是我用于测试的样本。我使用的是配备 i7-8565U CPU、Win 10 …

我正在尝试对以下函数进行矢量化以删除 sapply 循环。我正在计算累积偏度。

cskewness <- function(.x) {
  skewness <- function(.x) {
    sqrt(length(.x)) * sum((.x - mean(.x))^3) / (sum((.x - mean(.x))^2)^(3 / 2))
  }
  sapply(seq_along(.x), function(k, z) skewness(z[1:k]), z = .x)
}

我的代数没搞对。有这个是错误的：

skewness2 <- function(.x) {
  n <- length(.x)
  csum <- cumsum(.x)
  cmu <- csum / 1:length(.x)
  num <- cumsum(.x - cmu)^3
  den <- cumsum((.x - cmu)^2)^(3/2)
  sqrt(n) * num / den
}

正确的代码会产生：

x <- c(1,2,4,5,8)

> cskewness(x)
[1]       NaN 0.0000000 0.3818018 0.0000000 0.4082483
> skewness2(x)
[1]      NaN 1.000000 …

假设在MATLAB中,我有一个矩阵A,其元素为0或1.

如何以更快的矢量化方式获取每列的最后一个非零元素的索引向量？

我可以

[B, I] = max(cumsum(A));

和使用I,但有更快的方法吗？(我假设cumsum会花费一点时间甚至加0和1).

编辑:我想我矢量甚至比我需要快速的更多-福兹先生"循环是伟大的,但在MATLAB每个循环似乎花费了我很多的调试即使是快的时间.

我有一大段代码,其中一部分包含这段代码:

result = (nx * m_Lx + ny * m_Ly + m_Lz) / sqrt(nx * nx + ny * ny + 1);

我已经矢量化如下(一切都已经是float):

__m128 r = _mm_mul_ps(_mm_set_ps(ny, nx, ny, nx),
                      _mm_set_ps(ny, nx, m_Ly, m_Lx));
__declspec(align(16)) int asInt[4] = {
    _mm_extract_ps(r,0), _mm_extract_ps(r,1),
    _mm_extract_ps(r,2), _mm_extract_ps(r,3)
};
float (&res)[4] = reinterpret_cast<float (&)[4]>(asInt);
result = (res[0] + res[1] + m_Lz) / sqrt(res[2] + res[3] + 1);

结果是正确的; 但是,我的基准测试显示矢量化版本较慢:

• 非矢量化版本需要3750毫秒
• 矢量化版本需要4050毫秒
• 设置result到0直接(和完全除去这部分代码)降低了整个过程至2500毫秒

鉴于矢量版本只包含一个集SSE乘法(而不是四个单独的FPU乘法)的,为什么会慢？FPU确实比SSE快,或者这里有混淆变量吗？

(我在移动Core i5上.)

我需要找到一个向量的元素,这些元素少于它之后的一个或多个元素.在循环中很容易做到:

x = some_vector_values;
for m = 1 : length(x)
  if( any( x(m+1:end) > x(m) )
    do_such_and_such;
  end
end

但速度正在扼杀我.我正试图想出一个有效的解决办法但是我的空白.数组长度可以是数千种,我需要为许多不同的数组执行此操作.

我有两个具有形状(N,2,2)的3d阵列A和B,我想根据N轴在每个2x2矩阵上使用矩阵乘积.使用循环实现,它看起来像

C[i] = dot(A[i], B[i])

有没有办法在不使用循环的情况下做到这一点？我已经研究过tensordot,但还是无法让它工作.我想我可能会想要一些类似于tensordot(a, b, axes=([1,2], [2,1]))但却给我一个NxN矩阵的东西.

使用的numpy的可以减去的形状(3)阵列的优良广播规则v从一个形状(5,3)排列X以

X - v

结果是一个shape(5,3)数组,其中每一行i都是差异X[i] - v.

有没有办法减去形状(n,3)数组w,X以便从w整个数组中减去每一行X而不显式使用循环？

我有两个三维数组,a和b,并且想要找到b的2D子阵列,其中a的元素沿着第三轴具有最小值,即

a=n.random.rand(20).reshape((5,2,2))
b=n.arange(20).reshape((5,2,2))
c=n.argmin(a,2) #indices with minimal value of a
d=n.zeros_like(c) #the array I want
for i in range(5):
  for j in range(2):
     d[i,j] = b[i,j,c[i,j]] 

有没有办法可以在没有双循环的情况下获得这些值？

我知道这个答案: 在numpy数组中将min值替换为另一个 但是如果我想让它适用于我的3D数组,我必须做很多重塑操作 - 我想知道是否有更简单的东西.

我有以下numpy随机2D数组:

np.random.rand(30, 20)

我想迭代数组中的每个网格单元格.如果网格单元的值> 0.6,那么我想将多余的分配给其直接的8个相邻单元(在角网格单元的情况下,相邻单元的数量将更少).

应根据2个用户选择的规则之一重新分配超出部分:

1. 在8个邻居之间
2. 与每个邻居中的值成比例,即具有较高值的​​邻居变得更高

有没有办法在numpy不诉诸for循环的情况下完成此操作？