Covariance matrices estimation speed

[lkorczowski]

Hi,

I wanted to know if there is a specific reason why, in asr_calibrate, the following loop was used for the covariance estimation

line 158

% calculate the sample covariance matrices U (averaged in blocks of blocksize successive samples)
U = zeros(length(1:blocksize:S),C*C);
for k=1:blocksize
    range = min(S,k:blocksize:(S+k-1));
    U = U + reshape(bsxfun(@times,reshape(X(range,:),[],1,C),reshape(X(range,:),[],C,1)),size(U));
end

A faster and less memory intensive implementation would be something like this:

U2 = zeros(length(1:blocksize:S),C*C);
indx=1:blocksize:S;
for i=k:length(1:blocksize:S)
    indices=indx(k):min(S,indx(k)+blocksize-1);
    U2(k,:) = reshape(X(indices,:)' * X(indices,:),[1 C*C]);
end

Also, I see that that the default blocksize is 10 which seems very low for a correct empirical covariance estimation. Shouldn't we prefer a blocksize at least greater than the number of electrodes (or even blocksize>2*C) ?

Thanks!

[lkorczowski]

example for the two implementation:

S = 100001

C = 64

blocksize = 100

rASR implementation (U)
Elapsed time is 1.240739 seconds.
proposed implementation (U2)
Elapsed time is 0.077896 seconds.

copy/past the code to test yourself

S=100001
C=64
blocksize=100
X=randn(S,C);
% calculate the sample covariance matrices U (averaged in blocks of blocksize successive samples)
U = zeros(length(1:blocksize:S),C*C);
tic
for k=1:blocksize
    range = min(S,k:blocksize:(S+k-1));
    U = U + reshape(bsxfun(@times,reshape(X(range,:),[],1,C),reshape(X(range,:),[],C,1)),size(U));
end
toc

U2 = zeros(length(1:blocksize:S),C*C);

tic
indx=1:blocksize:S;
for k=1:length(1:blocksize:S)
    indices=indx(k):min(S,indx(k)+blocksize-1);
    U2(k,:) = reshape(X(indices,:)' * X(indices,:),[1 C*C]);
end
toc
sum(sum(abs(U-U2)))

comparison between the two implementation for different blocksize (I did every 10)

[s4rify]

Thank you very much for the proposed improvement!
To answer your initial question: there is no reason (for me) to compute the covariance matrix like in the original code, I just did not exchange the original code so far. I am checking out your suggestion and will let you know about the progress!