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, така че не е изключително голямо, но това изчисление се извършва на всяка времева стъпка на интегратор Runge Kutta. Освен това трябва да осредня статистическа съвкупност от произволни матрици 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(:)

Бенчмаркинг

Този раздел сравнява времената на изпълнение за двата подхода, изброени в тази публикация, първия тестван подход в публикацията на Shai и оригиналния подход, изброен във въпроса .

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

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. Относно бонуса, това би било страхотна идея, но в рутината Runge Kutta 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