Как применить функции corr2 в многомерных массивах в Matlab?

Допустим, у меня есть две матрицы A и B.

A = rand(4,5,3);
B = rand(4,5,6)

Я хочу применить функцию «corr2» для расчета коэффициентов корреляции.

corr2(A(:,:,1),B(:,:,1))
corr2(A(:,:,1),B(:,:,2))
corr2(A(:,:,1),B(:,:,3))
...
corr2(A(:,:,1),B(:,:,6))
...
corr2(A(:,:,2),B(:,:,1))
corr2(A(:,:,2),B(:,:,2))
...
corr2(A(:,:,3),B(:,:,6))

Как избежать использования циклов для создания такой векторизации?


person Chong Tang    schedule 23.10.2014    source источник


Ответы (2)


Взломан m-файл для corr2, чтобы создать кастомизированную векторизованную версию для работы с 3D-массивами. Здесь предлагаются два подхода с bsxfun (конечно же!)

Подход №1

szA = size(A);
szB = size(B);

a1 = bsxfun(@minus,A,mean(mean(A)));
b1 = bsxfun(@minus,B,mean(mean(B)));

sa1 = sum(sum(a1.*a1));
sb1 = sum(sum(b1.*b1));

v1 = reshape(b1,[],szB(3)).'*reshape(a1,[],szA(3));
v2 = sqrt(sb1(:)*sa1(:).');

corr3_out = v1./v2; %// desired output

corr3_out сохраняет результаты corr2 между всеми 3D-срезами A и B.

Таким образом, для A = rand(4,5,3), B = rand(4,5,6) у нас будет corr3_out как массив 6x3.

Подход № 2

Немного другой подход, чтобы сэкономить на нескольких вызовах sum и mean, используя вместо этого reshape -

szA = size(A);
szB = size(B);
dim12 = szA(1)*szA(2);

a1 = bsxfun(@minus,A,mean(reshape(A,dim12,1,[])));
b1 = bsxfun(@minus,B,mean(reshape(B,dim12,1,[])));

v1 = reshape(b1,[],szB(3)).'*reshape(a1,[],szA(3));
v2 = sqrt(sum(reshape(b1.*b1,dim12,[])).'*sum(reshape(a1.*a1,dim12,[])));

corr3_out = v1./v2; %// desired output

Бенчмаркинг

Код эталона -

%// Create random input arrays
N = 55; %// datasize scaling factor
A = rand(4*N,5*N,3*N);
B = rand(4*N,5*N,6*N);

%// Warm up tic/toc
for k = 1:50000
    tic(); elapsed = toc(); 
end

%// Run vectorized and loopy approach codes on the input arrays

%// 1. Vectorized approach
%//... solution code (Approach #2) posted earlier
%// clear variables used

%// 2. Loopy approach
tic
s_A=size(A,3);
s_B=size(B,3);
out1 = zeros(s_B,s_A);
for ii=1:s_A
    for jj=1:s_B
        out1(jj,ii)=corr2(A(:,:,ii),B(:,:,jj));
    end
end
toc

Полученные результаты -

-------------------------- With BSXFUN vectorized solution 
Elapsed time is 1.231230 seconds.
-------------------------- With loopy approach
Elapsed time is 139.934719 seconds.

Любители MATLAB-JIT проявляют любовь здесь! :)

person Divakar    schedule 23.10.2014
comment
Ага. это делает трюк. вы можете получить улучшение в 1000 раз для больших матриц - person ASantosRibeiro; 23.10.2014
comment
Если кому-то интересно, вот среда выполнения для циклического подхода, когда corr2 код встроен для того же размера данных, что и при тестировании - Elapsed time is 83.948572 seconds. Это позволяет избежать всех ненужных проверок ошибок и прочего. Я добавлю их в следующей редакции. - person Divakar; 23.10.2014
comment
Не могу дождаться, чтобы увидеть вашу следующую версию - person Chong Tang; 23.10.2014
comment
@ user2502865 О! Я имел в виду следующую редакцию (при необходимости, например, если вам нужны еще какие-либо изменения или правки). Но это готово для проверки! :) Вы пробовали коды? :) - person Divakar; 23.10.2014
comment
Я проверю этот код на огромных данных на следующей неделе. После этого я выберу лучший ответ на этот вопрос. - person Chong Tang; 25.10.2014
comment
@ user2502865 Теперь вы знаете размер этих огромных данных? - person Divakar; 26.10.2014
comment
Результатом может быть массив 30 000 * 300. - person Chong Tang; 26.10.2014
comment
@user2502865 user2502865 А каковы размеры A и B? Есть информация об этом? - person Divakar; 26.10.2014
comment
Пока еще не могу этого сказать. - person Chong Tang; 26.10.2014

Несколько примеров, но нет ничего лучше, чем циклы. Как говорит Дивакар в комментарии ниже, это не векторизованное решение.

КОД:

A = rand(4,5,1000);
B = rand(4,5,200);
s_A=size(A,3);
s_B=size(B,3);

%%% option 1
tic
corr_AB=cell2mat(arrayfun(@(indx1) arrayfun(@(indx2) corr2(A(:,:,indx1),B(:,:,indx2)),1:s_B),1:s_A,'UniformOutput',false));
toc

%%% option 2
tic
indx1=repmat(1:s_A,s_B,1);
indx1=indx1(:);
indx2=repmat(1:s_B,1,s_A);
indx2=indx2(:);
indx=[indx1,indx2];

corr_AB=arrayfun(@(i) corr2(A(:,:,indx(i,1)),B(:,:,indx(i,2))),1:size(indx,1));
toc

%%% option 3
tic
a=1;
for i=1:s_A
    for j=1:s_B
        corr_AB(a)=corr2(A(:,:,i),B(:,:,j));
        a=a+1;
    end
end
toc

ВЫХОД:

Elapsed time is 9.655696 seconds.
Elapsed time is 9.398979 seconds.
Elapsed time is 8.489744 seconds.
person ASantosRibeiro    schedule 23.10.2014
comment
Я проверяю это прямо сейчас. Любое объяснение вашего кода? - person Chong Tang; 23.10.2014
comment
некоторое редактирование, так как другое не работало должным образом. Я не вижу простого способа избежать петель для этого случая. для получения дополнительной информации читайте arrayfun - person ASantosRibeiro; 23.10.2014
comment
Большое спасибо. @АсантосРибейро - person Chong Tang; 23.10.2014
comment
Кажется, что ни один из циклов не лучше, чем циклы for, следует установленной тенденции - matlab-why">arrayfun может быть значительно медленнее, чем явный цикл в Matlab. Почему? По той же причине araryfun на самом деле не векторизованное решение. - person Divakar; 23.10.2014
comment
Это становится действительно фантастическим вопросом. ха-ха - person Chong Tang; 23.10.2014