Умножение 3D-матрицы на 2D-матрицу
https://stackoverflow.com/questions/1745299
-
20-09-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
Предположим, у меня есть АксВхС матрица X и БхД матрица Y.
Есть ли метод без цикла, с помощью которого я могу умножить каждое из С АксБ матрицы с Y?
Решение
Вы можете сделать это в одной строке, используя функции NUM2CELL сломать матрицу X в массив ячеек и СЕЛФАН для работы в ячейках:
Z = cellfun(@(x) x*Y,num2cell(X,[1 2]),'UniformOutput',false);
Результат Z это 1 по C массив ячеек, где каждая ячейка содержит A-by-D матрица.Если ты хочешь Z быть A-by-D-by-C матрицу, вы можете использовать КОТ функция:
Z = cat(3,Z{:});
ПРИМЕЧАНИЕ: Мое старое решение использовалось MAT2CELL вместо NUM2CELL, что было не так лаконично:
[A,B,C] = size(X);
Z = cellfun(@(x) x*Y,mat2cell(X,A,B,ones(1,C)),'UniformOutput',false);
Другие советы
Лично я предпочитаю, чтобы мой код был максимально кратким и читабельным.
Вот что я бы сделал, хотя это не соответствует вашему требованию «без петель»:
for m = 1:C
Z(:,:,m) = X(:,:,m)*Y;
end
Это приводит к А х Д х С матрица З.
И, конечно же, вы всегда можете заранее выделить Z, чтобы ускорить работу, используя Z = zeros(A,D,C);.
Вот однострочное решение (две, если вы хотите разделить на третье измерение):
A = 2;
B = 3;
C = 4;
D = 5;
X = rand(A,B,C);
Y = rand(B,D);
%# calculate result in one big matrix
Z = reshape(reshape(permute(X, [2 1 3]), [A B*C]), [B A*C])' * Y;
%'# split into third dimension
Z = permute(reshape(Z',[D A C]),[2 1 3]);
Следовательно, теперь: Z(:,:,i) содержит результат X(:,:,i) * Y
Объяснение:
Вышеизложенное может показаться запутанным, но идея проста.Сначала я начну с третьего измерения X и выполните вертикальную конкатенацию вдоль первого дима:
XX = cat(1, X(:,:,1), X(:,:,2), ..., X(:,:,C))
...трудность заключалась в том, что C является переменной, поэтому вы не можете обобщить это выражение, используя кот или верткат.Далее мы умножаем это на Y:
ZZ = XX * Y;
Наконец я разделил его обратно на третье измерение:
Z(:,:,1) = ZZ(1:2, :);
Z(:,:,2) = ZZ(3:4, :);
Z(:,:,3) = ZZ(5:6, :);
Z(:,:,4) = ZZ(7:8, :);
Итак, вы можете видеть, что для этого требуется только одно умножение матрицы, но вам нужно изменить форму Матрица до и после.
Я подхожу к той же проблеме, пытаясь найти наиболее эффективный метод.Я вижу примерно три подхода, за исключением использования внешних библиотек (т.е. mtimesx):
- Перебирать фрагменты 3D-матрицы
- волшебство повторения и перестановки
- умножение CellFun
Недавно я сравнил все три метода, чтобы определить, какой из них самый быстрый.Моя интуиция подсказывала, что победителем станет (2).Вот код:
% generate data
A = 20;
B = 30;
C = 40;
D = 50;
X = rand(A,B,C);
Y = rand(B,D);
% ------ Approach 1: Loop (via @Zaid)
tic
Z1 = zeros(A,D,C);
for m = 1:C
Z1(:,:,m) = X(:,:,m)*Y;
end
toc
% ------ Approach 2: Reshape+Permute (via @Amro)
tic
Z2 = reshape(reshape(permute(X, [2 1 3]), [A B*C]), [B A*C])' * Y;
Z2 = permute(reshape(Z2',[D A C]),[2 1 3]);
toc
% ------ Approach 3: cellfun (via @gnovice)
tic
Z3 = cellfun(@(x) x*Y,num2cell(X,[1 2]),'UniformOutput',false);
Z3 = cat(3,Z3{:});
toc
Все три подхода дали одинаковый результат (ух!), но, что удивительно, цикл оказался самым быстрым:
Elapsed time is 0.000418 seconds.
Elapsed time is 0.000887 seconds.
Elapsed time is 0.001841 seconds.
Обратите внимание, что время может сильно различаться от одного испытания к другому, и иногда (2) оказывается самым медленным.Эти различия становятся более существенными при увеличении объема данных.Но с много большие данные, (3) бьют (2).Метод цикла по-прежнему является лучшим.
% pretty big data...
A = 200;
B = 300;
C = 400;
D = 500;
Elapsed time is 0.373831 seconds.
Elapsed time is 0.638041 seconds.
Elapsed time is 0.724581 seconds.
% even bigger....
A = 200;
B = 200;
C = 400;
D = 5000;
Elapsed time is 4.314076 seconds.
Elapsed time is 11.553289 seconds.
Elapsed time is 5.233725 seconds.
Но метод цикла может быть медленнее, чем (2), если зацикленное измерение намного больше остальных.
A = 2;
B = 3;
C = 400000;
D = 5;
Elapsed time is 0.780933 seconds.
Elapsed time is 0.073189 seconds.
Elapsed time is 2.590697 seconds.
Таким образом, (2) в этом (возможно, крайнем) случае выигрывает с большим преимуществом.Возможно, не существует подхода, оптимального во всех случаях, но цикл по-прежнему довольно хорош и во многих случаях является лучшим.Это также лучший вариант с точки зрения читабельности.Отправляйтесь в путь!
Неа.Есть несколько способов, но всегда получается петля, прямая или косвенная.
Просто чтобы удовлетворить мое любопытство, зачем вам это вообще нужно?
Чтобы ответить на вопрос, и для удобства чтения см.:
- ндмульт, автор ajuanpi (Хуан Пабло Карбахал), 2013, GNU GPL
Вход
- 2 массива
- тусклый
Пример
nT = 100;
t = 2*pi*linspace (0,1,nT)’;
# 2 experiments measuring 3 signals at nT timestamps
signals = zeros(nT,3,2);
signals(:,:,1) = [sin(2*t) cos(2*t) sin(4*t).^2];
signals(:,:,2) = [sin(2*t+pi/4) cos(2*t+pi/4) sin(4*t+pi/6).^2];
sT(:,:,1) = signals(:,:,1)’;
sT(:,:,2) = signals(:,:,2)’;
G = ndmult (signals,sT,[1 2]);
Источник
Оригинальный источник.Я добавил встроенные комментарии.
function M = ndmult (A,B,dim)
dA = dim(1);
dB = dim(2);
# reshape A into 2d
sA = size (A);
nA = length (sA);
perA = [1:(dA-1) (dA+1):(nA-1) nA dA](1:nA);
Ap = permute (A, perA);
Ap = reshape (Ap, prod (sA(perA(1:end-1))), sA(perA(end)));
# reshape B into 2d
sB = size (B);
nB = length (sB);
perB = [dB 1:(dB-1) (dB+1):(nB-1) nB](1:nB);
Bp = permute (B, perB);
Bp = reshape (Bp, sB(perB(1)), prod (sB(perB(2:end))));
# multiply
M = Ap * Bp;
# reshape back to original format
s = [sA(perA(1:end-1)) sB(perB(2:end))];
M = squeeze (reshape (M, s));
endfunction
Я настоятельно рекомендую вам использовать Набор инструментов MMX из матлаба.Он может перемножать n-мерные матрицы максимально быстро.
Преимущества ММХ являются:
- Это легкий использовать.
- Умножить n-мерные матрицы (на самом деле он может умножать массивы двумерных матриц)
- Он выполняет другие матричные операции (транспонирование, квадратичное умножение, разложение Chol и многое другое)
- Оно использует Компилятор Си и многопотоковый расчет для ускорения.
Для решения этой проблемы вам просто нужно написать эту команду:
C=mmx('mul',X,Y);
вот эталон для всех возможных методов.Для более подробной информации обратитесь к этому вопрос.
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 <===================
Я бы подумал о рекурсии, но это единственный метод без цикла, который вы можете использовать.
Вы можете «развернуть» цикл, т.е. последовательно выписать все умножения, которые будут происходить в цикле.
Я хотел бы поделиться своим ответом на проблемы:
1) составление тензорного произведения двух тензоров (любой валентности);
2) осуществление сжатия двух тензоров по любому измерению.
Вот мои подпрограммы для первой и второй задач:
1) тензорное произведение:
function [C] = tensor(A,B)
C = squeeze( reshape( repmat(A(:), 1, numel(B)).*B(:).' , [size(A),size(B)] ) );
end
2) сокращение:Здесь A и B — тензоры, которые необходимо сжать по измерениям i и j соответственно.Длины этих размеров, конечно, должны быть равны.Для этого нет никакой проверки (это запутает код), но в остальном все работает хорошо.
function [C] = tensorcontraction(A,B, i,j)
sa = size(A);
La = length(sa);
ia = 1:La;
ia(i) = [];
ia = [ia i];
sb = size(B);
Lb = length(sb);
ib = 1:Lb;
ib(j) = [];
ib = [j ib];
% making the i-th dimension the last in A
A1 = permute(A, ia);
% making the j-th dimension the first in B
B1 = permute(B, ib);
% making both A and B 2D-matrices to make use of the
% matrix multiplication along the second dimension of A
% and the first dimension of B
A2 = reshape(A1, [],sa(i));
B2 = reshape(B1, sb(j),[]);
% here's the implicit implication that sa(i) == sb(j),
% otherwise - crash
C2 = A2*B2;
% back to the original shape with the exception
% of dimensions along which we've just contracted
sa(i) = [];
sb(j) = [];
C = squeeze( reshape( C2, [sa,sb] ) );
end
Есть критики?
