Быстрый и эффективный способ создания матрицы из серии продуктов

Question

Ax, Ay, Az: [N-by-N]

B = AA (диадическое произведение)

Это означает:

B(i,j)= [Ax(i,j);Ay(i,j);Az(i,j)]*[Ax(i,j) Ay(i,j) Az(i,j)]

B (i, j): матрица 3x3.Один из способов построить B:

N=2;
Ax=rand(N); Ay=rand(N); Az=rand(N);     %# [N-by-N]
t=1;
F=zeros(3,3,N^2);
for  i=1:N
    for j=1:N
        F(:,:,t)= [Ax(i,j);Ay(i,j);Az(i,j)]*[Ax(i,j) Ay(i,j) Az(i,j)];
        t=t+1;                          %# t is just a counter
    end
end
%# then we can write
B = mat2cell(F,3,3,ones(N^2,1));
B = reshape(B,N,N)'; 
B = cell2mat(B);

Есть ли более быстрый способ для больших N *.

Редактировать:

Спасибо за ваш ответ.(Это быстрее) Давайте положим: N = 2;Ax = [1 2; 3 4];Ay = [5 6; 7 8];Az = [9 10; 11 12];

B =

 1     5     9     4    12    20
 5    25    45    12    36    60
 9    45    81    20    60   100
 9    21    33    16    32    48
21    49    77    32    64    96
33    77   121    48    96   144

Выполнение:
???Ошибка при использовании ==> mtimes Размеры внутренней матрицы должны совпадать.

Если я напишу: P = Ai*Aj;, то

B  =

 7    19    31    15    43    71
23    67   111    31    91   151
39   115   191    47   139   231
10    22    34    22    50    78
34    78   122    46   106   166
58   134   210    70   162   254

То, что не зависит сверху A (:,:, 1) зависит от[Топор (1,1) Ай (1,1) Аз (1,1)]

Редактировать:

N=100;
Me      :Elapsed time is 1.614244 seconds.
gnovice :Elapsed time is 0.056575 seconds.
N=200;
Me      :Elapsed time is 6.044628 seconds.
gnovice :Elapsed time is 0.182455 seconds.
N=400;
Me      :Elapsed time is 23.775540 seconds.
gnovice :Elapsed time is 0.756682 seconds.
Fast!

rwong: B was not the same.

Редактировать:

После некоторой модификации для моегоприменение: по кодам gnovice

 1st code :  19.303310 seconds
 2nd code:  23.128920  seconds
 3rd code:  13.363585  seconds

Кажется, что любая функция, вызывающая, например, ceil, ind2sub ... делает thw циклы медленными и избегают, если это возможно.

symIndex было интересно!Спасибо.

Answer 1

Вот довольно простая и общая реализация, которая использует один цикл for для выполнения линейного индексирования и избегает работы с трехмерными переменными или изменения формы:

%# General solution:
%# ----------------
B = cell(N);
for index = 1:N^2
  A = [Ax(index) Ay(index) Az(index)];
  B{index} = A(:)*A;
end
B = cell2mat(B);

РЕДАКТИРОВАТЬ # 1:

В ответ на дополнительный вопрос о том, как уменьшить количество вычислений, выполняемых для симметричной матрицы B, вы можете использовать следующую модифицированную версию вышеуказанного кода:

%# Symmetric solution #1:
%# ---------------------
B = cell(N);
for index = find(tril(ones(N))).'  %'# Loop over the lower triangular part of B
  A = [Ax(index) Ay(index) Az(index)];
  B{index} = A(:)*A;
  symIndex = N*rem(index-1,N)+ceil(index/N);  %# Find the linear index for the
                                              %#   symmetric element
  if symIndex ~= index       %# If we're not on the main diagonal
    B{symIndex} = B{index};  %#   then copy the symmetric element
  end
end
B = cell2mat(B);

Однако, в таком случае вы можете получить более высокую производительность (или, по крайней мере, более простой код), отказавшись от линейного индексирования и используя два цикла for с индексированным индексированием, например, так:

%# Symmetric solution #2:
%# ---------------------
B = cell(N);
for c = 1:N    %# Loop over the columns
  for r = c:N  %# Loop over a subset of the rows
    A = [Ax(r,c) Ay(r,c) Az(r,c)];
    B{r,c} = A(:)*A;
    if r ~= c           %# If we're not on the main diagonal
      B{c,r} = B{r,c};  %#   then copy the symmetric element
    end
  end
end
B = cell2mat(B);

РЕДАКТИРОВАНИЕ № 2:

Второе симметричное решение может быть еще более ускорено путем перемещения диагонального вычисления за пределы внутреннего цикла (устраняя необходимость в условном выражении) и перезаписи A с результатом A(:)*A, чтобы мы могли обновить симметричную ячейку введите B{c,r}, используя A вместо B{r,c} (т.е. обновление одной ячейки содержимым другой, по-видимому, несет некоторые дополнительные издержки):

%# Symmetric solution #3:
%# ---------------------
B = cell(N);
for c = 1:N    %# Loop over the columns
  A = [Ax(c,c) Ay(c,c) Az(c,c)];
  B{c,c} = A(:)*A;
  for r = c+1:N  %# Loop over a subset of the rows
    A = [Ax(r,c) Ay(r,c) Az(r,c)];
    A = A(:)*A;
    B{r,c} = A;
    B{c,r} = A;
  end
end
B = cell2mat(B);

А вот некоторые временные результаты для 3 симметричных решений, использующих следующие примерные симметричные матрицы Ax, Ay и Az:

N = 400;
Ax = tril(rand(N));     %# Lower triangular matrix
Ax = Ax+triu(Ax.',1);  %'# Add transpose to fill upper triangle
Ay = tril(rand(N));     %# Lower triangular matrix
Ay = Ay+triu(Ay.',1);  %'# Add transpose to fill upper triangle
Az = tril(rand(N));     %# Lower triangular matrix
Az = Az+triu(Az.',1);  %'# Add transpose to fill upper triangle

%# Timing results:
%# --------------
%# Solution #1 = 0.779415 seconds
%# Solution #2 = 0.704830 seconds
%# Solution #3 = 0.325920 seconds

Большое ускорение для решения 3 обусловлено главным образом обновлением содержимого ячейки B локальной переменной A вместо обновления одной ячейки содержимым другой. Другими словами, копирование содержимого ячейки с помощью таких операций, как B{c,r} = B{r,c};, несет больше затрат, чем я ожидал.

Answer 2

A = cat(3, Ax, Ay, Az);   % [N-by-N-by-3]
F = zeros(3, 3, N^2);
for i = 1:3, 
  for j = 1:3,
    Ai = A(:,:,i);
    Aj = A(:,:,j);
    P = Ai(:) .* Aj(:);
    F(i,j,:) = reshape(P, [1, 1, N^2]);
  end
end

%# then we can write
B = mat2cell(F,3,3,ones(N^2,1));
B = reshape(B,N,N)'; 
B = cell2mat(B);