эффективно вычислять 3D матрицу внешних продуктов - MATLAB

Question

Предположим, у меня есть матрица элементов следующим образом:

A = reshape(1:25, 5, 5) 

A = 

1  6 11 16 21 
2  7 12 17 22 
3  8 13 18 23 
4  9 14 19 24 
5 10 15 20 25 

Я хотел бы эффективно вычислять 3D матрицу внешних продуктов, так что я ^й срезом этой выходной матрицы является внешнее произведение столбца i ^thA с самим собой. Внешний продукт между двумя векторами u и v равен u*v.', если u и v являются обоими векторами-столбцами.

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

B(:,:,1) = A(:,1) * A(:,1).'; 
B(:,:,2) = A(:,2) * A(:,2).'; 
     ... 
     ... 
B(:,:,5) = A(:,5) * A(:,5).'; 

Мой текущий метод состоит в следующем. Я пытался делать это таким образом, используя arrayfun и cell2mat:

cellmatr = arrayfun(@(x) A(:,x) * A(:,x).', 1:size(A,2), 'uni', 0); 
out = reshape(cell2mat(cellmatr), size(A,1), size(A,1), size(A,2)); 

Я просто перебираем линейный массив индекса между 1 и столько же столбцов, которые мы имеем в A, и для каждого элемента в этом массиве, я получить доступ к соответствующей колонке и вычислить внешний продукт. Таким образом, выход будет давать 1D сетку ячеек, которые затем преобразуются обратно в 2D-матрицу, а затем преобразуются в трехмерную матрицу, чтобы найти 3D-матрицу внешних продуктов.

Однако для больших матриц это довольно медленно. Я также попытался заменить матричный продукт kron (т. Е. kron(A(:,x), A(:,x))) внутри моего вызова arrayfun, но это все еще довольно медленно для моих целей.

---

Кто-нибудь знает об эффективном способе вычисления этой 3D-матрицы внешних продуктов таким образом?

Answer 1

Это лишь незначительное улучшение по сравнению с Divakar's answer. Это немного быстрее, потому что он заменяет переставлять 3D-массив с перестановочны 2D-массива:

B = bsxfun(@times, permute(A, [1 3 2]), permute(A, [3 1 2])); 

Answer 2

Это -

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

Answer 3

Чтобы констатировать очевидное, вы пробовали просто для цикла:

[m,n] = size(A); 
B = zeros(m,m,n); 
for i=1:n 
    B(:,:,i) = A(:,i) * A(:,i).'; 
end 

Вы будете удивлены, насколько быстро конкурентно это.