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 × k矩阵,β是k × 1随机向量.我们通过蒙特卡罗整合计算这种期望-让小号平的测试,计算EXP(X*试用版)对每个平局,然后取平均值.通常我们得到m> k> j的数据,我们使用非常大的s.在这个例子中,我简单地让X成为1的矩阵.

我使用了3个版本的矢量化bsxfun(),它们的区别在于X和β的形状:

矢量化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总是最快？

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

编辑:这篇文章的早期版本不正确.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)

x我们的最终目标是使用和betausing中的“消除”k matrix-multiplication。因此，我们可以将其“推”到最后并k重新整形为保持行的二维，即 ( j * m , k )，然后与( k , s ) 执行矩阵乘法，得到 ( j * m , s )。然后可以将乘积重新整形为 3D 数组 (j, m, s) 并执行元素指数运算，结果为。xpermutekbetaexp_xBeta

现在，如果最终目标是y，即沿 的第三维求平均值exp_xBeta，则相当于计算沿矩阵乘积 (j * m, s ) 的行的平均值，然后重塑为 ( j , m ) 直接联系我们y。