Correlation coefficients between two matrices to find intercorrelation

Question

I am trying to calculate Pearson coefficients between all pair combinations of my variables of all my samples.

Say i have an m*n matrix where m are the variables and n are the samples

i want to calculate for each variable of my data what is the correlation to every other variable.

So, i managed to do that with nested loops:

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

and it works, but takes too much time to run

I am trying to find a way to run the corrcoef function on a row basis, meaning maybe to create an additional matrix with repmat of the base values and correlate to the X data using some FUN function.

Could not figure out how to use the fun with inputs from to arrays, running between individuals lines/columns

help will be appreciated

@GameOfThrows doesn't corrcoef only return a 2x2 matrix in this case? your suggestion wouldn't produce the correct number of numerical outputs — hiandbaii
@hiandbaii yes you are right, that was a silly comment, I have deleted it. I think arrayfun is the way to solve this — GameOfThrows

Divakar Divakar · Accepted Answer · 2016-01-11T20:04:25

This post involves a bit of hacking, so bear with it!

Stage #0 To start off, we have -

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 From the documentation of corrcoef in its source code :

If C is the covariance matrix, C = COV(X), then CORRCOEF(X) is the matrix whose (i,j)'th element is : C(i,j)/SQRT(C(i,i)*C(j,j)).

After hacking into the code of covariance, we see that for the default case of one input, the covariance formula is simply -

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

Thus, mixing the two definitions and putting them into the problem at hand, we have -

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

Stage #2 This is the final stage where we use the real fun aka bsxfun to kill all loops, like so -

%// 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.'));

Benchmarking

This section compares the original approach against the proposed one and also verifies the output. Here's the benchmarking code -

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(:)))

Output -

---------- 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

Correlation coefficients between two matrices to find intercorrelation

1 Answers

Benchmarking