Matlab bsxfun() - 什么解释了沿不同维度扩展时的性能差异？

Question

在我的工作（计量经济学/统计学）中，我经常需要将不同大小的矩阵相乘，然后对生成的矩阵执行额外的操作。我一直依赖bsxfun() 对代码进行矢量化处理，总的来说我发现它比repmat() 更有效。但我不明白的是，为什么有时bsxfun() 的性能在沿不同维度扩展矩阵时会非常不同。

考虑这个具体的例子：

x      = ones(j, k, m);
beta   = rand(k, m, s);

exp_xBeta   = zeros(j, m, s);
for im = 1 : m
    for is = 1 : s
        xBeta                = x(:, :, im) * beta(:, im, is);
        exp_xBeta(:, im, is) = exp(xBeta);
    end
end

y = mean(exp_xBeta, 3);

上下文：

我们有来自 m 个市场的数据，在每个市场中，我们想要计算 exp(X * beta) 的期望值，其中 X 是j x k 矩阵，beta 是 k x 1 随机向量.我们通过 monte-carlo 积分来计算这个期望值 - 对 beta 进行 s 次绘制，为每次绘制计算 exp(X * beta)，然后取平均值。通常我们使用 m > k > j 获取数据，并且我们使用非常大的 s。在这个例子中，我只是让 X 成为一个矩阵。

我使用bsxfun() 进行了 3 个版本的矢量化，它们的不同之处在于 X 和 beta 的形状：

矢量化 1

x1      = x;                                    % size [ j k m 1 ]
beta1   = permute(beta, [4 1 2 3]);             % size [ 1 k m s ]

tic
xBeta       = bsxfun(@times, x1, beta1);
exp_xBeta   = exp(sum(xBeta, 2));
y1          = permute(mean(exp_xBeta, 4), [1 3 2 4]);   % size [ j m ]
time1       = toc;

矢量化 2

x2      = permute(x, [4 1 2 3]);                % size [ 1 j k m ]
beta2   = permute(beta, [3 4 1 2]);             % size [ s 1 k m ]

tic
xBeta       = bsxfun(@times, x2, beta2);
exp_xBeta   = exp(sum(xBeta, 3));
y2          = permute(mean(exp_xBeta, 1), [2 4 1 3]);   % size [ j m ]
time2       = toc;

矢量化 3

x3      = permute(x, [2 1 3 4]);                % size [ k j m 1 ]
beta3   = permute(beta, [1 4 2 3]);             % size [ k 1 m s ]

tic
xBeta       = bsxfun(@times, x3, beta3);
exp_xBeta   = exp(sum(xBeta, 1));
y3          = permute(mean(exp_xBeta, 4), [2 3 1 4]);    % size [ j m ]
time3       = toc;

这就是他们的表现（通常我们使用 m > k > j 获取数据，并且我们使用了非常大的 s）：

j = 5，k = 15，m = 100，s = 2000：

For-loop version took 0.7286 seconds.
Vectorized version 1 took 0.0735 seconds.
Vectorized version 2 took 0.0369 seconds.
Vectorized version 3 took 0.0503 seconds.

j = 10，k = 15，m = 150，s = 5000：

For-loop version took 2.7815 seconds.
Vectorized version 1 took 0.3565 seconds.
Vectorized version 2 took 0.2657 seconds.
Vectorized version 3 took 0.3433 seconds.

j = 15，k = 35，m = 150，s = 5000：

For-loop version took 3.4881 seconds.
Vectorized version 1 took 1.0687 seconds.
Vectorized version 2 took 0.8465 seconds.
Vectorized version 3 took 0.9414 seconds.

为什么版本 2 始终是最快的版本？最初，我认为性能优势是因为 s 设置为维度 1，Matlab 可能能够更快地计算，因为它以列优先顺序存储数据。但 Matlab 的分析器告诉我，计算该平均值所花费的时间相当微不足道，并且在所有 3 个版本中或多或少相同。 Matlab 大部分时间都在评估 bsxfun() 行，这也是 3 个版本中运行时差异最大的地方。

有没有想过为什么版本 1 总是最慢而版本 2 总是最快？

我在这里更新了我的测试代码： Code

编辑：这篇文章的早期版本不正确。 beta 的大小应为 (k, m, s)。

Answer 1

bsxfun 当然是向量化事物的好工具之一，但如果你能以某种方式引入matrix-multiplication，那将是最好的方法，如matrix multiplications are really fast on MATLAB。

看来在这里你可以像这样使用matrix-multiplication 来获取exp_xBeta -

[m1,n1,r1] = size(x);
n2 = size(beta,2);
exp_xBeta_matmult = reshape(exp(reshape(permute(x,[1 3 2]),[],n1)*beta),m1,r1,n2)

或者直接获取y如下图-

y_matmult = reshape(mean(exp(reshape(permute(x,[1 3 2]),[],n1)*beta),2),m1,r1)

说明

为了更详细地解释它，我们的尺寸为 -

x      : (j, k, m)
beta   : (k, s)

我们的最终目标是使用matrix-multiplication 来“消除”来自x 和beta 的k。因此，我们可以将x 中的k“推”到permute 的末尾，并重塑为保持k 作为行的二维，即（j * m，k），然后执行矩阵乘法beta ( k , s ) 给我们 ( j * m , s )。然后可以将产品重新整形为 3D 数组（ j ， m ， s ）并执行元素指数，这将是 exp_xBeta。

现在，如果最终目标是y，即沿exp_xBeta 的第三维求均值，则相当于沿矩阵乘积 (j * m, s ) 然后再整形为 ( j , m ) 直接得到我们y。

Answer 2

今天早上我又做了一些实验。这似乎与Matlab以列主顺序存储数据的事实有关。

在做这些实验时，我还添加了矢量化版本 4，它做同样的事情，但订购的维度与版本 1-3 略有不同。

---

回顾一下，以下是所有 4 个版本中 x 和 beta 的排序方式：

矢量化 1：

x       :   (j, k, m, 1)
beta    :   (1, k, m, s)

矢量化 2：

x       :   (1, j, k, m)
beta    :   (s, 1, k, m)

矢量化 3：

x       :   (k, j, m, 1)
beta    :   (k, 1, m, s)

矢量化 4：

x       :   (1, k, j, m)
beta    :   (s, k, 1, m)

代码：bsxfun_test.m

---

这段代码中开销最大的两个操作是：

(a)xBeta = bsxfun(@times, x, beta);

(b)exp_xBeta = exp(sum(xBeta, dimK));

其中dimK 是k 的维度。

在 (a) 中，bsxfun() 必须沿 s 的维度扩展 x 和沿 j 的维度扩展 beta。当s 远大于其他维度时，我们应该看到向量化 2 和 4 的一些性能优势，因为它们将s 指定为第一个维度。

j = 100; k = 100; m = 100; s = 1000;

Vectorized version 1 took 2.4719 seconds.
Vectorized version 2 took 2.1419 seconds.
Vectorized version 3 took 2.5071 seconds.
Vectorized version 4 took 2.0825 seconds.

如果s 微不足道而k 很大，那么矢量化3 应该是最快的，因为它将k 放入维度1：

j = 10; k = 10000; m = 100; s = 1;

Vectorized version 1 took 0.0329 seconds.
Vectorized version 2 took 0.1442 seconds.
Vectorized version 3 took 0.0253 seconds.
Vectorized version 4 took 0.1415 seconds.

如果我们在最后一个示例中交换 k 和 j 的值，则向量化 1 将变得最快，因为 j 被分配到维度 1：

j = 10000; k = 10; m = 100; s = 1;

Vectorized version 1 took 0.0316 seconds.
Vectorized version 2 took 0.1402 seconds.
Vectorized version 3 took 0.0385 seconds.
Vectorized version 4 took 0.1608 seconds.

但一般来说，当k 和j 接近时，j > k 并不一定意味着向量化 1 比向量化 3 快，因为 (a) 和 (b) 中执行的操作不同。

在实践中，我经常需要使用s >>>> m > k > j 运行计算。在这种情况下，似乎在矢量化 2 或 4 中对它们进行排序会得到最好的结果：

    j = 10; k = 30; m = 100; s = 5000;

Vectorized version 1 took 0.4621 seconds.
Vectorized version 2 took 0.3373 seconds.
Vectorized version 3 took 0.3713 seconds.
Vectorized version 4 took 0.3533 seconds.

    j = 15; k = 50; m = 150; s = 5000;

Vectorized version 1 took 1.5416 seconds.
Vectorized version 2 took 1.2143 seconds.
Vectorized version 3 took 1.2842 seconds.
Vectorized version 4 took 1.2684 seconds.

---

要点：如果bsxfun() 必须沿着比其他维度大得多的维度扩展，则将该维度分配给维度 1！

Answer 3

请参考其他question 和answer

如果您要使用bsxfun 处理不同维度的矩阵，请确保矩阵的最大维度保留在第一维度中。

这是我的小示例测试：

%// Inputs
%// Taking one very big and one small vector, so that the difference could be seen clearly
a = rand(1000000,1);
b = rand(1,5);

%//---------------- testing with inbuilt function
%// preferred orientation [1]
t1 = timeit(@() bsxfun(@times, a, b))

%// not preferred [2]
t2 = timeit(@() bsxfun(@times, b.', a.'))

%//---------------- testing with anonymous function
%// preferred orientation [1]
t3 = timeit(@() bsxfun(@(x,y) x*y, a, b))

%// not preferred [2]
t4 = timeit(@() bsxfun(@(x,y) x*y, b.', a.'))

^{[1] 首选方向      -    较大的维度作为第一维度
[2] 不喜欢                 -    较小的维度作为第一维度}

小注：这四种方法给出的输出都是一样的，尽管它们的维度可能不同。

结果：

t1 =
0.0461

t2 =
0.0491

t3 =
0.0740

t4 =
7.5249

>> t4/t3
ans =
101.6878

---

Method 3 大约比Method 4 快100 倍

---

总结： 尽管内置函数的首选方向和不喜欢的方向之间的差异很小， 对于匿名函数，差异变得巨大。因此，最好使用更大的维度作为维度 1。