向量化：矩阵数组逐个元素相乘答案

【问题标题】：Vectorization: matrix array multiplication element wise one by one向量化：矩阵数组逐个元素相乘
【发布时间】：2017-07-09 02:53:02
【问题描述】：

我有一个矩阵：

R = [0 -1;1 0];

array = 1:1:10;

还有x0 = [2;1]

如何在没有循环的情况下以最有效的方式获取另一个数组？

array2 = [expm(1*R) expm(2*R) expm(3*R) .... expm(10*R)];

那我想获取 array3 的维度为 2 x 10，这样：

array3 = [expm(1*R)*x0 expm(2*R)*x0 expm(3*R)*x0 .... expm(10*R)*x0];

【问题讨论】：

R 是向量还是矩阵？这是没有意义的，您的 array3 可以是 2 x 20 矩阵或 2 x 10 x 2 矩阵。请澄清。
R 是一个矩阵。 array3 是 2 乘 10。
正如我所说，这没有任何意义。 array 具有维度 1x10，R 具有维度 2x2 生成的矩阵 array2 或 array3 不能少于 2x2x10 = 40 元素。请给出一个带有数字的确切示例输出。
@thewaywewalk 关于array2的大小，expm(1*R)的大小是2x2。在array2 中水平连接了 10 个这样的矩阵，使其大小为 2x20。而关于array3的大小，expm(1*R)*x0的大小是2x1，在array3中有10个这样的矩阵水平连接，使其大小为2x10。
@SardarUsama 所以expm(1*R)*x0 应该是真正的矩阵乘法而不是元素？关于不清楚的问题的标题。

标签： matlab vectorization

【解决方案1】：

来自wikipedia：

如果一个矩阵是对角矩阵，它的指数可以通过对主对角线上的每个元素取幂来获得。

假设可以从{1*R, 2*R,...} 创建块对角矩阵，则可以获取其指数并将其重新整形为[2 * n]，并且可以乘以x0。但是它的性能可能比 for 循环差。

R = [0 -1;1 0];
array = 1:1:10;
x0 = [2;1]
n = numel(array);
result = reshape(expm(kron(spdiags(array.',0,n,n),R))*repmat(x0,n,1),2,[]);

对于小尺寸（少于 70 个元素）的array，全矩阵更高效：

result = reshape(expm(kron(diag(array),R))*repmat(x0,n,1),2,[]);

【讨论】：

您应该尝试以某种方式将kron 替换为bsxfun，我认为这是您解决方案的瓶颈。（见基准）。不过想法不错！ +1

【解决方案2】：

好吧，我看到您拥有的矩阵 R 是 2x2。如果始终是 2x2，则可以使用以下函数 (Wikipedia) 来计算指数：

function output = expm2d(A)
% Assuming t = 1 from Evaluation by Laurent series (https://en.wikipedia.org/wiki/Matrix_exponential#Evaluation_by_Laurent_series)
s = trace(A) / 2;
q = sqrt(-det(A - s*eye(size(A))));
output = exp(s) * ((cosh(q) - s * sinh(q) / q) * eye(size(A)) + (sinh(q) * A / q));
end

使用thewaywewalk提供的优秀比较功能，得到如下结果：

使用expm时：

>> bench
ans =
    0.0181 %// rahnema
    0.1075 %// thewaywewalk arrayfun
    0.1139 %// thewaywewalk accumarray

使用expm2d时：

>> bench
ans =
    0.0048 %// rahnema
    0.0161 %// thewaywewalk arrayfun
    0.0222 %// thewaywewalk accumarray

如您所见，对 2d 矩阵使用该函数会导致运行时间减少 10 倍。当然，当 R 不是 2x2 时，这个就不能用了。

编辑：将expm2d 用于A = 1:100 时：

>> bench
ans =
    0.1379 %// rahnema
    0.1415 %// thewaywewalk arrayfun
    0.1756 %// thewaywewalk accumarray

【讨论】：

【解决方案3】：

我仍然不知道你的问题是否正确。以下是两个未完全矢量化但相当快的解决方案：

R = [0 -1;1 0];
A = 1:1:10;
x0 = [2;1];

%// option 1
temp = arrayfun(@(x) (expm(R*x)*x0).', A, 'uni', 0);
array3 = vertcat( temp{:} )

%// option 2
temp = accumarray( (1:numel(A)).', A(:), [], @(x) {(expm(R*x)*x0).'})
array3 = vertcat( temp{:} )

基准测试

我没有考虑Leander's Answer，因为它不计算array3：

function [t] = bench()
    R = [0 -1;1 0];
    A = 1:1:10;
    x0 = [2;1];

    % functions to compare
    fcns = {
        @() compare1(A,R,x0);
        @() compare2(A,R,x0);
        @() compare3(A,R,x0);
    };

    % timeit
    t = zeros(3,1);
    for ii = 1:100;
        t = t + cellfun(@timeit, fcns);
    end
end

function array3 = compare1(A,R,x0)  %rahnema1
    n = numel(A);
    array3 = reshape(expm(kron(diag(A),R))*repmat(x0,n,1),2,[])
end
function array3 = compare2(A,R,x0)  %thewaywewalk 1
    temp = arrayfun(@(x) (expm(R*x)*x0).', A, 'uni', 0);
    array3 = vertcat( temp{:} )
end
function array3 = compare3(A,R,x0)  %thewaywewalk 2
    temp = accumarray( (1:numel(A)).', A(:), [], @(x) {(expm(R*x)*x0).'});
    array3 = vertcat( temp{:} )
end

结果：

为A = 1:1:10;

0.1006   %// rahnema
0.2831   %// thewaywewalk arrayfun
0.3103   %// thewaywewalk accumarray

由于kron 对于大型数组来说变得非常慢，A = 1:1:100; 的基准测试结果发生了变化：

4.0068   %// rahnema
1.8045   %// thewaywewalk arrayfun
2.4257   %// thewaywewalk accumarray

【讨论】：