logo标签列表

矩阵乘法中的bsxfun实现

performancematlabmatrixmatrix-multiplicationbsxfun
7
一如既往地希望从您那里学到更多，我希望能够得到以下代码的帮助。
我需要完成以下任务：
1）我有一个向量：
x = [1 2 3 4 5 6 7 8 9 10 11 12]

2) 以及一个矩阵：
A =[11    14    1
    5     8    18
    10    8    19
    13    20   16]

我需要能够将x中的每个值与A中的每个值相乘，这意味着：
new_matrix = [1* A
              2* A
              3* A
               ...
              12* A]

这将给我一个大小为（12 * m x n）的new_matrix，假设A（mxn）。在这种情况下，（12 * 4x3） 我该如何使用Matlab中的bsxfun来实现这一点？这种方法比for-loop更快吗？
关于我的for-loop，我需要一些帮助......当循环运行时，我无法存储每个"new_matrix" :(
for i=x
new_matrix = A.*x(i)
end

提前感谢！

编辑：在解决方案提供后

第一个解决方案

clear all
clc
x=1:0.1:50;
A = rand(1000,1000);
tic
val = bsxfun(@times,A,permute(x,[3 1 2]));
out = reshape(permute(val,[1 3 2]),size(val,1)*size(val,3),[]);
toc

输出：

Elapsed time is 7.597939 seconds.

第二个解决方案
clear all
clc
x=1:0.1:50;
A = rand(1000,1000);
tic
Ps = kron(x.',A);
toc

输出：

Elapsed time is 48.445417 seconds.
for循环可以通过预定义大小为(12*m,n)的new_matrix来实现，就像你自己说的那样，然后使用索引告诉new_matrix你想要保存这些元素的位置，例如在你上面给出的代码中，new_matrix(((i-1)*12+1):(i*12))) = A.*x(i)，我在这里写了它，所以不确定它是否有效。 - The Minion
谢谢@Minion，我会检查它是否有效并告诉你！ - Sergio Haram
@Minion 它基本上是有效的，但我得到了一些介于 1*new_matrix、2*new_matrix、3*new_matrix 等等之间的东西，还有其他一些计算，我无法确定它们来自哪里。 - Sergio Haram
1
@SergioHaram 谢谢您发布这个问题！希望对那些对 bsxfun 感兴趣的人有所帮助。 - Divakar
2
酷！一些基准测试结果！！感谢您发布这些！ - Divakar
5个回答
17

将x发送到第三维，这样在使用bsxfun与A相乘时，单例展开将生效，扩展乘积结果至第三维。然后进行bsxfun乘法运算 -

val = bsxfun(@times,A,permute(x,[3 1 2]))

现在，val是一个3D矩阵，期望的输出是通过第三维沿列拼接的2D矩阵。下面实现了这个目标 -

out = reshape(permute(val,[1 3 2]),size(val,1)*size(val,3),[])

希望这讲得清楚！请传播bsxfun的知识！woo!! :)

1
这应该运行得很好。至少对我的测试来说是这样的。但我真的不理解这段代码。为什么要排列？你介意编辑一下你的答案并解释一下吗？ - The Minion
@Divakar 好的，我理解了将 x 扩展到第三维，这是因为 A 的 dim 是这样的，对吗？在第二行中，您正在将 mapping 返回到所需的 dim。但是，如果我将数字 [3 1 2] 更改为任何其他顺序，代码是否仍然有效？你的第二行中有 [1 3 2](!)？关于第二行：size(val,3) 和 [] 是什么意思？ - Sergio Haram
谢谢大家！我担心我发了一个非常无聊的问题...(也许是吧)，但是这个小时我学到了很多！！！ - Sergio Haram
1
@Divakar 很好地使用了 bsxfun。 - Luis Mendo
@Divakar 抱歉我没有表达清楚。如果我的矩阵 A 是 (1000x1000)，而我的向量 x 是 (500x1)（只是举个例子），如果我在 val = bsxfun(@times,A,permute(x,[3 1 2])) 中只有三个变量 (3 1 2)，那么单例将会超出 dim，这样说对吗？ - Sergio Haram
显示剩余11条评论
10

kron 函数正好可以实现这一点：

kron(x.',A)
Kronecker product（克罗内克积）！不错 @Luis Mendo！ - Sergio Haram
1
谢谢大家！我担心我发了一个非常无聊的问题...(也许是吧)，但是这个小时我学到了很多！！！ - Sergio Haram
1
@Divakar 是的，很棒的工具。如果你需要做相同的事情，比如求和而不是乘积，可以看一下 kron 的代码：它使用 meshgrid 生成所有组合，没有更多的东西。 - Luis Mendo
3

以下是我对目前提到的方法以及我自己添加的一些方法的基准测试：

function [t,v] = testMatMult()
    % data
    %{
    x = [1 2 3 4 5 6 7 8 9 10 11 12];
    A = [11 14 1; 5 8 18; 10 8 19; 13 20 16];
    %}
    x = 1:50;
    A = randi(100, [1000,1000]);

    % functions to test
    fcns = {
        @() func1_repmat(A,x)
        @() func2_bsxfun_3rd_dim(A,x)
        @() func2_forloop_3rd_dim(A,x)
        @() func3_kron(A,x)
        @() func4_forloop_matrix(A,x)
        @() func5_forloop_cell(A,x)
        @() func6_arrayfun(A,x)
    };

    % timeit
    t = cellfun(@timeit, fcns, 'UniformOutput',true);

    % check results
    v = cellfun(@feval, fcns, 'UniformOutput',false);
    isequal(v{:})
    %for i=2:numel(v), assert(norm(v{1}-v{2}) < 1e-9), end
end

% Amro
function B = func1_repmat(A,x)
    B = repmat(x, size(A,1), 1);
    B = bsxfun(@times, B(:), repmat(A,numel(x),1));
end

% Divakar
function B = func2_bsxfun_3rd_dim(A,x)
    B = bsxfun(@times, A, permute(x, [3 1 2]));
    B = reshape(permute(B, [1 3 2]), [], size(A,2));
end

% Vissenbot
function B = func2_forloop_3rd_dim(A,x)
    B = zeros([size(A) numel(x)], 'like',A);
    for i=1:numel(x)
        B(:,:,i) = x(i) .* A;
    end
    B = reshape(permute(B, [1 3 2]), [], size(A,2));
end

% Luis Mendo
function B = func3_kron(A,x)
    B = kron(x(:), A);
end

% SergioHaram & TheMinion
function B = func4_forloop_matrix(A,x)
    [m,n] = size(A);
    p = numel(x);
    B = zeros(m*p,n, 'like',A);
    for i=1:numel(x)
        B((i-1)*m+1:i*m,:) = x(i) .* A;
    end
end

% Amro
function B = func5_forloop_cell(A,x)
    B = cell(numel(x),1);
    for i=1:numel(x)
        B{i} = x(i) .* A;
    end
    B = cell2mat(B);
    %B = vertcat(B{:});
end

% Amro
function B = func6_arrayfun(A,x)
    B = cell2mat(arrayfun(@(xx) xx.*A, x(:), 'UniformOutput',false));
end

我的电脑上的结果：

>> t
t =
    0.1650    %# repmat (Amro)
    0.2915    %# bsxfun in the 3rd dimension (Divakar)
    0.4200    %# for-loop in the 3rd dim (Vissenbot)
    0.1284    %# kron (Luis Mendo)
    0.2997    %# for-loop with indexing (SergioHaram & TheMinion)
    0.5160    %# for-loop with cell array (Amro)
    0.4854    %# arrayfun (Amro)

（这些时间可能会在不同的运行中略有变化，但这应该让我们了解这些方法的比较）

请注意，对于较大的输入，其中一些方法将导致内存不足错误（例如基于 repmat 的解决方案很容易耗尽内存）。其他方法针对较大尺寸会显着变慢，但不会因为内存耗尽而出错（例如 kron 解决方案）。

我认为，在这种情况下， bsxfun 方法 func2_bsxfun_3rd_dim 或直接的for-loop func4_forloop_matrix （感谢MATLAB JIT）是最好的解决方案。

当然，您可以更改上述基准参数（ x 和 A 的大小），并得出自己的结论 :)

2

补充一种可选方案，你可以使用cellfun来实现你想要的功能。以下是一个示例（略有修改）：

x = randi(2, 5, 3)-1;
a = randi(3,3);
%// bsxfun 3D (As implemented in the accepted solution)
val = bsxfun(@and, a, permute(x', [3 1 2])); %//'
out = reshape(permute(val,[1 3 2]),size(val,1)*size(val,3),[]);
%// cellfun (My solution)
val2 = cellfun(@(z) bsxfun(@and, a, z), num2cell(x, 2), 'UniformOutput', false);
out2 = cell2mat(val2); % or use cat(3, val2{:}) to get a 3D matrix equivalent to val and then permute/reshape like for out
%// compare
disp(nnz(out ~= out2));

两者都会得到完全相同的结果。

如果想了解更多使用cellfun的技巧和信息，请参考：http://matlabgeeks.com/tips-tutorials/computation-using-cellfun/

还有这个网址：https://dev59.com/3XI-5IYBdhLWcg3wpqIK#1746422

0

如果您的向量x长度为12，矩阵大小为3x4，我认为在时间方面使用其中之一不会有太大的变化。如果您正在处理更大的矩阵和向量，则可能会成为问题。

因此，首先我们想要将向量与矩阵相乘。在for循环方法中，这将给出类似以下的内容：

s = size(A);
new_matrix(s(1),s(2),numel(x)) = zeros;   %This is for pre-allocating. If you have a big vector or matrix, this will help a lot time efficiently.

for i = 1:numel(x)
    new_matrix(:,:,i)= A.*x(i)
end

这将为您提供3D矩阵，其中每个第三维度都是您乘法的结果。如果这不是您要寻找的内容，我将添加另一种解决方案，可能更适用于更大的矩阵和向量，并且更加高效。

1
非常感谢！我确实得到了我需要的矩阵，但是我只需要一个包含使用您的代码时获得的所有结果的 new_matrix。是的，我正在处理非常大的矩阵，在这里，我只是给出了一个简单的例子来让其他人理解我的问题。 - Sergio Haram
是的，当我看到Divakar的答案后我注意到了！我看错了，以为它是一个12xNxM矩阵，而不是12*NxM！Divakar的答案对我来说看起来很不错！ - Vissenbot
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接