Каков быстрый способ вычисления корреляции столбца по столбцу в Matlab

У меня есть две очень большие матрицы (60x25000), и я хотел бы вычислить корреляцию между столбцами только между двумя матрицами. Например:

corrVal(1) = corr(mat1(:,1), mat2(:,1);
corrVal(2) = corr(mat1(:,2), mat2(:,2);
...
corrVal(i) = corr(mat1(:,i), mat2(:,i);

Для меньших матриц я могу просто использовать:

   colCorr = diag( corr( mat1, mat2 ) );

но это не работает для очень больших матриц, так как у меня заканчивается память. Я подумал о том, чтобы разделить матрицы на части для вычисления корреляций, а затем объединить результаты, но вычисление корреляции между комбинациями столбцов кажется бесполезным, что мне на самом деле неинтересно.

Есть ли быстрый способ напрямую вычислить то, что мне интересно?

Изменить. Раньше я использовал цикл, но его просто замедлили:

mat1 = rand(60,5000);
mat2 = rand(60,5000);
nCol = size(mat1,2);
corrVal = zeros(nCol,1);

tic;
for i = 1:nCol
    corrVal(i) = corr(mat1(:,i), mat2(:,i));
end
toc;

Это занимает ~ 1 секунду

tic;
corrVal = diag(corr(mat1,mat2));
toc;

Это занимает ~ 0,2 секунды.


matlab correlation matrix
person slayton    schedule 13.02.2012    source источник
comment
Я отредактировал ваш пост; проверьте, правильно ли это.   -  person Jacob    schedule 13.02.2012
comment
Кроме того, что не так с очевидным циклом for?   -  person Jacob    schedule 13.02.2012
comment
редактирование правильное, спасибо! Также цикл замедляет   -  person slayton    schedule 13.02.2012
comment
Я сделал еще одно изменение. Кроме того, на моем ПК цикл занял ~ 1,7 с, а версия diag все еще работает (более минуты).   -  person Jacob    schedule 13.02.2012
comment
Хорошо, я уменьшил матрицу до 60x500, а версии loop и diag заняли ~ 0,17 с и ~ 16,7 с соответственно.   -  person Jacob    schedule 13.02.2012
comment
хм ... вот интересно, какая версия Matlab у вас работает и сколько у вас ядер? Я использую 2011b, и у меня 8 ядер, которые могут давать преимущество на моей машине. Когда выполняется диагностика, я вижу всплеск активности на большинстве своих ядер.   -  person slayton    schedule 13.02.2012

Ответы (2)


arrow_upward
16
arrow_downward

Я могу получить увеличение скорости x100, вычислив его вручную.

An=bsxfun(@minus,A,mean(A,1)); %%% zero-mean
Bn=bsxfun(@minus,B,mean(B,1)); %%% zero-mean
An=bsxfun(@times,An,1./sqrt(sum(An.^2,1))); %% L2-normalization
Bn=bsxfun(@times,Bn,1./sqrt(sum(Bn.^2,1))); %% L2-normalization
C=sum(An.*Bn,1); %% correlation

Вы можете сравнить, используя этот код:

A=rand(60,25000);
B=rand(60,25000);

tic;
C=zeros(1,size(A,2));
for i = 1:size(A,2)
    C(i)=corr(A(:,i), B(:,i));
end
toc; 

tic
An=bsxfun(@minus,A,mean(A,1));
Bn=bsxfun(@minus,B,mean(B,1));
An=bsxfun(@times,An,1./sqrt(sum(An.^2,1)));
Bn=bsxfun(@times,Bn,1./sqrt(sum(Bn.^2,1)));
C2=sum(An.*Bn,1);
toc
mean(abs(C-C2)) %% difference between methods

Вот время вычислений:

Elapsed time is 10.822766 seconds.
Elapsed time is 0.119731 seconds.

Разница между двумя результатами очень мала:

mean(abs(C-C2))

ans =
  3.0968e-17

РЕДАКТИРОВАТЬ: объяснение

bsxfun выполняет операцию "столбец за столбцом" (или "строка за строкой" в зависимости от ввода).

An=bsxfun(@minus,A,mean(A,1));

Эта строка удалит (@minus) среднее значение каждого столбца (mean(A,1)) для каждого столбца _9 _... Таким образом, в основном это делает столбцы A с нулевым средним.

An=bsxfun(@times,An,1./sqrt(sum(An.^2,1)));

В этой строке умножается (@ раз) каждый столбец на значение, обратное его норме. Таким образом, это делает их L-2 нормализованными.

После того, как столбцы имеют нулевое среднее значение и нормализованы по L2, для вычисления корреляции вам просто нужно произвести скалярное произведение каждого столбца An на каждый столбец B. Таким образом, вы умножаете их поэлементно An.*Bn, а затем суммируете каждый столбец: sum(An.*Bn);.

person Oli    schedule 13.02.2012
comment
вау, это быстро. Не могли бы вы дать мне быстрое объяснение, почему это работает? - person slayton; 13.02.2012
comment
Спасибо uuuuu за отличный ответ @Oli - person Christina; 12.04.2016

arrow_upward
1
arrow_downward

Я думаю, что очевидный цикл может быть достаточно хорош для вашего размера проблемы. На моем ноутбуке для выполнения следующих действий требуется менее 6 секунд:

A = rand(60,25000);
B = rand(60,25000);
n = size(A,1);
m = size(A,2);

corrVal = zeros(1,m);
for k=1:m
    corrVal(k) = corr(A(:,k),B(:,k));
end
person Ian Hincks    schedule 13.02.2012
comment
Подождите, что мне не хватает? diag (corr (A, B)); занимает у меня более 10 секунд. - person Ian Hincks; 13.02.2012
comment
Сможете ли вы запустить diag(corr(A,B)) на таких больших матрицах? Я получаю ошибку нехватки памяти, когда пытаюсь запустить эти матрицы - person slayton; 13.02.2012
comment
@slayton: Не думаю, что ты сможешь. Кроме того, кажется, что версия цикла быстрее! - person Jacob; 13.02.2012
comment
@slayton У меня 8 ГБ оперативной памяти, и она заполняет около 5 ГБ, но это работает. - person Ian Hincks; 13.02.2012