两个矩阵之间的相关系数以找到相关性

Question

我正在尝试计算所有样本的变量的所有对组合之间的 Pearson 系数。

假设我有一个 m*n 矩阵，其中 m 是变量，n 是样本

我想为我的数据中的每个变量计算与其他所有变量的相关性。

所以，我设法通过嵌套循环来做到这一点：

X = rand[1000 100];
for i = 1:1000
base = X(i, :);
    for j = 1:1000
    target = X(j, :);
    correlation = corrcoef(base, target);
    correlation = correlation(2, 1);
    corData(1, j) = correlation
    end
totalCor(i, :) = corData
end

它可以工作，但运行时间太长

我正在尝试找到一种方法来逐行运行 corrcoef 函数，这意味着可能会创建一个带有基值 repmat 的附加矩阵，并使用一些 FUN 函数与 X 数据相关联。

无法弄清楚如何使用来自数组的输入的乐趣，在个人行/列之间运行

我们将不胜感激

Answer 1

这篇文章涉及到一点黑客攻击，请耐心等待！

第 0 阶段首先，我们有 -

for i = 1:N
    base = X(i, :);
    for j = 1:N
        target = X(j, :);
        correlation = corrcoef(base, target);
        correlation = correlation(2, 1)
        corData(1, j) = correlation;
    end
end

Stage #1来自corrcoef源代码中的文档：

如果 C 是协方差矩阵，C = COV(X)，那么 CORRCOEF(X) 是 第 (i,j) 个元素为：C(i,j)/SQRT(C(i,i)*C(j,j)) 的矩阵。

在破解covariance的代码后，我们看到对于一个输入的默认情况，协方差公式很简单-

[m,n] = size(x);
xc = bsxfun(@minus,x,sum(x,1)/m);
xy = (xc' * xc) / (m-1);

因此，混合这两个定义并将它们放入手头的问题中，我们有 -

m = size(X,2);
for i = 1:N
    base = X(i, :);
    for j = 1:N
        target = X(j, :);
        BT = [base(:) target(:)];
        xc = bsxfun(@minus,BT,sum(BT,1)/m);
        C = (xc' * xc) / (m-1); %//'
        corData = C(2,1)/sqrt(C(2,2)*C(1,1))
    end
end

第 2 阶段这是我们使用 真正的乐趣 aka bsxfun 来杀死所有循环的最后阶段，就像这样 -

%// Broadcasted subtract of each row by the average of it.
%// This corresponds to "xc = bsxfun(@minus,BT,sum(BT,1)/m)"
p1 = bsxfun(@minus,X,mean(X,2));

%// Get pairs of rows from X and get the dot product. 
%// Thus, a total of "N x N" such products would be obtained.
p2 = sum(bsxfun(@times,permute(p1,[1 3 2]),permute(p1,[3 1 2])),3);

%// Scale them down by "size(X,2)-1". 
%// This was for the part : "C = (xc' * xc) / (m-1)".
p3 = p2/(size(X,2)-1);

%// "C(2,2)" and "C(1,1)" are diagonal elements from "p3", so store them.
dp3 = diag(p3);

%// Get "sqrt(C(2,2)*C(1,1))" by broadcasting elementwise multiplication 
%// of "dp3". Finally do elementwise division of "p3" by it.
totalCor_out = p3./sqrt(bsxfun(@times,dp3,dp3.'));

基准测试

本节将原始方法与提议的方法进行比较，并验证输出。这是基准测试代码 -

disp('---------- With original approach')
tic
X = rand(1000,100);
corData = zeros(1,1000);
totalCor = zeros(1000,1000);
for i = 1:1000
    base = X(i, :);
    for j = 1:1000
        target = X(j, :);
        correlation = corrcoef(base, target);
        correlation = correlation(2, 1);
        corData(1, j) = correlation;
    end
    totalCor(i, :) = corData;
end
toc

disp('---------- With the real fun aka BSXFUN')
tic
p1 = bsxfun(@minus,X,mean(X,2));
p2 = sum(bsxfun(@times,permute(p1,[1 3 2]),permute(p1,[3 1 2])),3);
p3 = p2/(size(X,2)-1);
dp3 = diag(p3);
totalCor_out = p3./sqrt(bsxfun(@times,dp3,dp3.')); %//'
toc

error_val = max(abs(totalCor(:)-totalCor_out(:)))

输出 -

---------- With original approach
Elapsed time is 186.501746 seconds.
---------- With the real fun aka BSXFUN
Elapsed time is 1.423448 seconds.
error_val =
    4.996e-16