Matlab: векторизация суммы 4D-матриц

Мне нужно выполнить в MATLAB следующий расчет:

где w и v — векторы с N элементами, а A — четырехмерная матрица (N^4 элементов). Это может быть достигнуто с помощью следующего педантичного кода:

N=10;
A=rand(N,N,N,N);
v=rand(N,1);
w=zeros(N,1);

for pp=1:N
  for ll=1:N
    for mm=1:N
      for nn=1:N
        w(pp)=w(pp)+A(pp,ll,mm,nn)*v(ll)*v(mm)*conj(v(nn));
      end
    end
  end
end

что очень медленно. Есть ли способ векторизовать такую ​​сумму в MATLAB?


matlab performance nested-loops matrix vectorization
person Juan Sebastian Totero    schedule 09.03.2015    source источник
comment
С какими размерами данных вы на самом деле имеете дело? Итак, каково ваше настоящее N?   -  person Divakar    schedule 09.03.2015
comment
в реальных расчетах N должно быть порядка ~100, поэтому не очень велико, но этот расчет выполняется на каждом временном шаге интегратора Рунге-Кутта. Кроме того, мне нужно усреднить статистическую совокупность случайных матриц A, поэтому векторизация могла бы оказать значительную помощь.   -  person Juan Sebastian Totero    schedule 09.03.2015
comment
А * в v* указывает на транспонирование?   -  person kkuilla    schedule 09.03.2015
comment
@kkuilla Conjugate Наверное.   -  person Divakar    schedule 09.03.2015
comment
является ли A разреженной матрицей?   -  person Shai    schedule 09.03.2015
comment
@kkuilla: да, это сложная операция сопряжения   -  person Juan Sebastian Totero    schedule 09.03.2015
comment
@Shai: это плотная матрица   -  person Juan Sebastian Totero    schedule 09.03.2015

Ответы (2)


arrow_upward
6
arrow_downward

Подход №1

С несколькими reshape и умножение матриц -

A1 = reshape(A,N^3,N)*conj(v)
A2 = reshape(A1,N^2,N)*v
w = reshape(A2,N,N)*v

Подход № 2

С одним bsxfun , reshape и matrix-multiplication -

A1 = reshape(A,N^3,N)*conj(v)
vm = bsxfun(@times,v,v.')
w = reshape(A1,N,N^2)*vm(:)

Бенчмаркинг

В этом разделе сравнивается время выполнения для двух подходов, перечисленных в этом сообщении, первый проверенный подход в сообщении Шай и исходный подход, указанный в вопросе. .

Код сравнительного анализа

N=100;
A=rand(N,N,N,N);
v=rand(N,1);

disp('----------------------------------- With Original Approach')
tic
%// .... Code from the original post   ...//
toc

disp('----------------------------------- With Shai Approach #1')
tic
s4 = sum( bsxfun( @times, A, permute( conj(v), [4 3 2 1] ) ), 4 ); 
s3 = sum( bsxfun( @times, s4, permute( v, [3 2 1] ) ), 3 );
w2 = s3*v; 
toc

disp('----------------------------------- With Divakar Approach #1')
tic
A1 = reshape(A,N^3,N)*conj(v);
A2 = reshape(A1,N^2,N)*v;
w3 = reshape(A2,N,N)*v;
toc

disp('----------------------------------- With Divakar Approach #2')
tic
A1 = reshape(A,N^3,N)*conj(v);
vm = bsxfun(@times,v,v.');
w4 = reshape(A1,N,N^2)*vm(:);
toc

Результаты выполнения

----------------------------------- With Original Approach
Elapsed time is 4.604767 seconds.
----------------------------------- With Shai Approach #1
Elapsed time is 0.334667 seconds.
----------------------------------- With Divakar Approach #1
Elapsed time is 0.071905 seconds.
----------------------------------- With Divakar Approach #2
Elapsed time is 0.058877 seconds.

Выводы

Второй подход в этом посте, похоже, дает примерно 80x ускорение по сравнению с исходным подходом.

person Divakar    schedule 09.03.2015
comment
Спасибо! Особенно для бенчмаркинга! Умный способ использовать bsxfun. Что касается бонуса, это было бы отличной идеей, но в программе Рунге-Кутта v является решением в момент времени n, а w — в момент времени n+1, поэтому я не могу его предварительно вычислить. :) - person Juan Sebastian Totero; 09.03.2015

arrow_upward
4
arrow_downward

Вы можете попробовать использовать bsxfun.

Предполагая, что v - это векторы-столбцы N на 1 (в противном случае перестановки должны быть немного изменены).

% sum over n (4th dim)
s4 = sum( bsxfun( @times, A, permute( conj(v), [4 3 2 1] ) ), 4 );

Теперь промежуточный результат всего N на N на N.

% sum over m (3rd dim)
s3 = sum( bsxfun( @times, s4, permute( v, [3 2 1] ) ), 3 )

Переходим к последней сумме

% sum over l (2nd dim)
w = s3*v;

Если подумать, не рассматривали ли вы возможность использования dot в его мультидим версия? Я не проверял, но должно работать (возможно небольшие исправления).

s4 = dot( A, permute( conj(v), [4 3 2 1] ), 4 );
s3 = dot( s4, permute( v, [3 2 1] ), 3 );
w = s3*v;

person Shai    schedule 09.03.2015
comment
Хорошее использование bsxfun должно быть здесь эффективным! - person Divakar; 09.03.2015
comment
Спасибо, ваш ответ был очень полезен! - person Juan Sebastian Totero; 09.03.2015