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

Вот довольно простая и общая реализация, в которой используется один цикл 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);

А вот некоторые результаты синхронизации для трех симметричных решений с использованием следующих образцов симметричных матриц 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};, влечет за собой больше накладных расходов, чем я ожидал.