МАТЛАБ: Как векторно-умножить два массива матриц? - PullRequest
Новая
       39

МАТЛАБ: Как векторно-умножить два массива матриц?

8 голосов
/ 05 июля 2011

У меня есть два 3-мерных массива, первые два измерения которых представляют матрицы, а последнее рассчитывает через пространство параметров, в качестве простого примера возьмем 

A = repmat([1,2; 3,4], [1 1 4]);

(но предположим, что A(:,:,j) отличаетсяза каждый j).Как можно легко выполнить умножение на матрицу за j двух таких матриц-массивов A и B? 

C = A; % pre-allocate, nan(size(A,1), size(B,2)) would be better but slower
for jj = 1:size(A, 3)
  C(:,:,jj) = A(:,:,jj) * B(:,:,jj);
end

, безусловно, выполняет свою работу, но если третье измерение больше похоже на 1e3элементы это очень медленно, так как он не использует векторизацию MATLAB.Так есть ли более быстрый способ?

Ответы [ 5 ]

6 голосов
/ 05 июля 2011

Я провел несколько временных тестов, самый быстрый способ для 2x2xN - это вычисление матричных элементов: 

C = A;
C(1,1,:) = A(1,1,:).*B(1,1,:) + A(1,2,:).*B(2,1,:);
C(1,2,:) = A(1,1,:).*B(1,2,:) + A(1,2,:).*B(2,2,:);
C(2,1,:) = A(2,1,:).*B(1,1,:) + A(2,2,:).*B(2,1,:);
C(2,2,:) = A(2,1,:).*B(1,2,:) + A(2,2,:).*B(2,2,:);

В общем случае выясняется, что цикл for на самом деле самый быстрый (не забудьте предварительно выделить C!).

Если у вас уже есть результат в виде матрицы ячеек, хотя, используя cellfun , это самый быстрый выбор, он также быстрее, чем зацикливание элементов ячейки: 

C = cellfun(@mtimes, A, B, 'UniformOutput', false);

Однако необходимость вызова num2cell first (Ac = num2cell(A, [1 2])) и cell2mat для случая 3d-массива тратит слишком много времени.

Вот некоторые моменты, которые я выбрал для случайного набора 2 x 2 x 1e4: 

 array-for: 0.057112
 arrayfun : 0.14206
 num2cell : 0.079468
 cell-for : 0.033173
 cellfun  : 0.025223
 cell2mat : 0.010213
 explicit : 0.0021338

Явное относится к использованию прямого вычисления матричных элементов 2 x 2, см. Ниже. Результат аналогичен для новых случайных массивов, cellfun является самым быстрым, если num2cell не требуется раньше и нет ограничений на 2x2xN. Для общих 3d-массивов циклы по третьему измерению - действительно самый быстрый выбор. Вот временной код: 

n = 2;
m = 2;
l = 1e4;

A = rand(n,m,l);
B = rand(m,n,l);

% naive for-loop:
tic
%Cf = nan(n,n,l);
Cf = A;
for jl = 1:l
    Cf(:,:,jl) = A(:,:,jl) * B(:,:,jl);
end;
disp([' array-for: ' num2str(toc)]);

% using arrayfun:
tic
Ca = arrayfun(@(k) A(:,:,k)*B(:,:,k), 1:size(A,3), 'UniformOutput',false);
Ca = cat(3,Ca{:});
disp([' arrayfun : ' num2str(toc)]);

tic
Ac = num2cell(A, [1 2]);
Bc = num2cell(B, [1 2]);
disp([' num2cell : ' num2str(toc)]);

% cell for-loop:
tic
Cfc = Ac;
for jl = 1:l
    Cfc{jl} = Ac{jl} * Bc{jl};
end;
disp([' cell-for : ' num2str(toc)]);

% using cellfun:
tic
Cc = cellfun(@mtimes, Ac, Bc, 'UniformOutput', false);
disp([' cellfun  : ' num2str(toc)]);

tic
Cc = cell2mat(Cc);
disp([' cell2mat : ' num2str(toc)]);

tic
Cm = A;
Cm(1,1,:) = A(1,1,:).*B(1,1,:) + A(1,2,:).*B(2,1,:);
Cm(1,2,:) = A(1,1,:).*B(1,2,:) + A(1,2,:).*B(2,2,:);
Cm(2,1,:) = A(2,1,:).*B(1,1,:) + A(2,2,:).*B(2,1,:);
Cm(2,2,:) = A(2,1,:).*B(1,2,:) + A(2,2,:).*B(2,2,:);
disp([' explicit : ' num2str(toc)]);

disp(' ');
4 голосов
/ 14 июля 2011

Вот мой контрольный тест, сравнивающий методы, упомянутые в @ TobiasKienzler ответе. Я использую функцию TIMEIT , чтобы получить более точные значения времени. 

function [t,v] = matrixMultTest()
    n = 2; m = 2; p = 1e5;
    A = rand(n,m,p);
    B = rand(m,n,p);

    %# time functions
    t = zeros(5,1);
    t(1) = timeit( @() func1(A,B,n,m,p) );
    t(2) = timeit( @() func2(A,B,n,m,p) );
    t(3) = timeit( @() func3(A,B,n,m,p) );
    t(4) = timeit( @() func4(A,B,n,m,p) );
    t(5) = timeit( @() func5(A,B,n,m,p) );

    %# check the results
    v = cell(5,1);
    v{1} = func1(A,B,n,m,p);
    v{2} = func2(A,B,n,m,p);
    v{3} = func3(A,B,n,m,p);
    v{4} = func4(A,B,n,m,p);
    v{5} = func5(A,B,n,m,p);
    assert( isequal(v{:}) )
end

%# simple FOR-loop
function C = func1(A,B,n,m,p)
    C = zeros(n,n,p);
    for k=1:p
        C(:,:,k) = A(:,:,k) * B(:,:,k);
    end
end

%# ARRAYFUN
function C = func2(A,B,n,m,p)
    C = arrayfun(@(k) A(:,:,k)*B(:,:,k), 1:p, 'UniformOutput',false);
    C = cat(3, C{:});
end

%# NUM2CELL/FOR-loop/CELL2MAT
function C = func3(A,B,n,m,p)
    Ac = num2cell(A, [1 2]);
    Bc = num2cell(B, [1 2]);
    C = cell(1,1,p);
    for k=1:p
        C{k} = Ac{k} * Bc{k};
    end;
    C = cell2mat(C);
end

%# NUM2CELL/CELLFUN/CELL2MAT
function C = func4(A,B,n,m,p)
    Ac = num2cell(A, [1 2]);
    Bc = num2cell(B, [1 2]);
    C = cellfun(@mtimes, Ac, Bc, 'UniformOutput', false);
    C = cell2mat(C);
end

%# Loop Unrolling
function C = func5(A,B,n,m,p)
    C = zeros(n,n,p);
    C(1,1,:) = A(1,1,:).*B(1,1,:) + A(1,2,:).*B(2,1,:);
    C(1,2,:) = A(1,1,:).*B(1,2,:) + A(1,2,:).*B(2,2,:);
    C(2,1,:) = A(2,1,:).*B(1,1,:) + A(2,2,:).*B(2,1,:);
    C(2,2,:) = A(2,1,:).*B(1,2,:) + A(2,2,:).*B(2,2,:);
end

Результаты: 

>> [t,v] = matrixMultTest();
>> t
t =
      0.63633      # FOR-loop
      1.5902       # ARRAYFUN
      1.1257       # NUM2CELL/FOR-loop/CELL2MAT
      1.0759       # NUM2CELL/CELLFUN/CELL2MAT
      0.05712      # Loop Unrolling

Как я объяснил в комментариях, простой цикл FOR является лучшим решением (за исключением разматывания цикла в последнем случае, что выполнимо только для этих небольших матриц размером 2 на 2).

3 голосов
/ 09 сентября 2015

Я настоятельно рекомендую вам использовать MMX toolbox из Matlab.Он может умножать n-мерные матрицы как можно быстрее.

Преимущества MMX :

  1. Это просто использовать.
  2. Умножение n-мерных матриц (фактически оно может умножать массивы двумерных матриц)
  3. Выполняет другие операции с матрицами (транспонировать,Квадратичное умножение, декомпозиция Chol и многое другое)
  4. Для ускорения используется C компилятор и многопоточность .

Для этогопроблема, вам просто нужно написать эту команду: 

C=mmx('mul',A,B);

Я добавил следующую функцию в ответ @ Amro 

%# mmx toolbox
function C=func6(A,B,n,m,p)
    C=mmx('mul',A,B);
end

Я получил этот результат для n=2,m=2,p=1e5: 

    1.6571 # FOR-loop
    4.3110 # ARRAYFUN
    3.3731 # NUM2CELL/FOR-loop/CELL2MAT
    2.9820 # NUM2CELL/CELLFUN/CELL2MAT
    0.0244 # Loop Unrolling
    0.0221 # MMX toolbox  <===================

Я использовал код @ Amro для запуска теста.

1 голос
/ 22 октября 2013

По моему опыту, еще более быстрый метод состоит в использовании точечного умножения и суммирования по трехмерной матрице.Следующая функция, функция z_matmultiply (A, B), умножает две трехмерные матрицы, которые имеют одинаковую глубину.Умножение точек выполняется максимально параллельным образом, поэтому вы можете проверить скорость этой функции и сравнить ее с другими при большом количестве повторений. 

function C = z_matmultiply(A,B)

[ma,na,oa] = size(A);
[mb,nb,ob] = size(B);

%preallocate the output as we will do a loop soon
C = zeros(ma,nb,oa);

%error message if the dimensions are not appropriate
if na ~= mb || oa ~= ob
    fprintf('\n z_matmultiply warning: Matrix Dimmensions Inconsistent \n')
else

% if statement minimizes for loops by looping the smallest matrix dimension 
if ma > nb
    for j = 1:nb
        Bp(j,:,:) = B(:,j,:);
        C(:,j,:) = sum(A.*repmat(Bp(j,:,:),[ma,1]),2);
    end
else
    for i = 1:ma
        Ap(:,i,:) = A(i,:,:);
        C(i,:,:) = sum(repmat(Ap(:,i,:),[1,nb]).*B,1);
    end 
end

end
1 голос
/ 05 июля 2011

Один из методов заключается в создании разреженной матрицы 2Nx2N и встраивании по диагонали матриц 2x2 как для A, так и для B. Сделайте произведение с разреженными матрицами и возьмите результат с немного хитрой индексацией и измените его до 2x2xN.

Но я сомневаюсь, что это будет быстрее, чем простой цикл.

...