Matlab中计算列与列之间相关性的快速方法是什么？

Question

我有两个非常大的矩阵（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秒钟的时间

- slayton

---

我对你的帖子进行了编辑，请检查是否正确。 - Jacob

还有，显然的for循环有什么问题吗？ - Jacob

编辑是正确的，谢谢！另外，循环速度太慢了。 - slayton

我又做了一个更改。此外，在我的电脑上，循环花费了约1.7秒，而“diag”版本仍在运行（已经超过一分钟）。 - Jacob

好的，我将矩阵缩小为60x500，循环和对角线版本分别需要约0.17秒和16.7秒。 - Jacob

嗯... 这很有趣，你运行的 MATLAB 版本是多少，有几个核心？我正在运行 2011b 版本，有 8 个核心，在我的机器上可能会给 diag 带来优势。当 diag 执行时，我看到大多数核心的活动都有所增加。 - slayton

Answer 1

我可以通过手工计算获得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));

这行代码将从A的每一列中减去该列的平均值 (mean(A,1))...所以基本上它使得A的每一列都是零均值。

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

这行代码通过每列的范数倒数进行乘法(@times)，使它们变为L-2归一化。

一旦列具有零均值和L2归一化，要计算相关性，只需将An的每一列与B的每一列进行点积。因此，您可以对它们进行逐元素乘法An.*Bn，然后对每一列进行求和：sum(An.*Bn);。

- Oli

---

哇，这真快。你能给我一个快速的解释为什么这个有效吗？ - slayton

非常感谢 @Oli 的出色回答！ - Christina

Answer 2

我认为对于你的问题规模来说，显然的循环可能已经足够好了。在我的笔记本电脑上，以下操作不到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

- Ian Hincks

---

等一下，我错过了什么？对我来说，diag(corr(A,B)); 花费了超过10秒钟的时间。 - Ian Hincks

你能在这么大的矩阵上运行 diag(corr(A,B)) 吗？当我尝试在这些矩阵上运行时，会出现内存错误。 - slayton

@slayton：我认为你不能这样做。此外，循环版本似乎更快！ - Jacob

@slayton 我有8GB的RAM，它占用了大约5GB的空间，但是它还能正常工作。 - Ian Hincks