在我的工作领域(计量经济学/统计学(中，我经常需要将不同大小的矩阵相乘，然后对生成的矩阵执行额外的操作。我一直依靠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 随机向量。我们通过蒙特卡罗积分来计算这个期望 - 进行 s 次抽奖，计算每次抽奖的 exp(X * beta(，然后取平均值。通常，我们获取 m> k> j 的数据，并且我们使用非常大的 s。在这个例子中，我只是让X成为1的矩阵。
我使用 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.
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.
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 总是最快的，有没有想过？
我在这里更新了我的测试代码：法典
编辑：这篇文章的早期版本不正确。 beta的大小应为 (k, m, s) .
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)
我们的最终目标是使用x中的"消除"k，并使用matrix-multiplication beta。因此，我们可以用permute将x中的k"推"到末尾，并重新塑造为行的 2D 保持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。
今天早上我又做了一些实验。这似乎与 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！
参考其他问题和答案
如果要使用 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。