Optimise/replace this for loop used in calculation of mutual information

So I previously posted on here with regards to some code given HERE. Thanks to the help of Jan Simons I have managed to eliminate one of the nested for loops to speed up the code pretty significantly. However, due to the dimensionality of the project and future projects I am working on I still need to optimise the code. I ran the profiler and the following block of code seems to be bottleneck:

Eps = zeros(nObs, 1);
Nn = zeros(nObs, 1);

nx1 = zeros(nObs, 1);
ny1 = zeros(nObs, 1);
nx2 = zeros(nObs, 1);
ny2 = zeros(nObs, 1);
for i = 1:nObs

    dxSample = dx(i, :);
    dxSample(i) = [];

    dySample = dy(i, :);
    dySample(i) = [];

    dzSample = dz(i, :);
    dzSample(i) = [];

    [EpsSample, NnSample] = sort(dzSample, 'ascend');

    Eps(i) = EpsSample(k);
    Nn(i) = NnSample(k);

    nx1(i) = sum(dxSample < Eps(i));
    ny1(i) = sum(dySample < Eps(i));

    nx2(i) = sum(dxSample <= Eps(i));
    ny2(i) = sum(dySample <= Eps(i));

 end

Note that the k is typically quite low: between 3-10, nObs is around 4364. Is there any way of speeding this up? Maybe even eliminating the for loop altogether if possible? Any help is very much appreciated!

[Greg]

What are dx, dy, and dz? Where does k come from?

Also, what release of MATLAB are you running? Implicit expansion is a life-saver.

[Jan]

@Greg: See https://www.mathworks.com/matlabcentral/answers/377090-speed-up-for-loop-in-this-code-for-calculating-mutual-information-maybe-using-gpu-computing#answer_300083. To my surprise replacing

 for i = j+1:n
   dy(i,j) = sqrt(sum(Y(:, i) - Yj ) .^ 2));
 end

by

 for i = j+1:n
   dy(i,j) = sqrt(sum(bsxfun(@minus, Y(:, i), Yj) .^ 2));
 end

gave a speedup of 10% in R2016b for the complete function. Implicit expansion is so much nicer, but not fast here.

Maybe these are valid test inputs (positive, symmetric, zero on the diagonal):

dx = rand(4364, 4364);
dy = rand(4364, 4364);
dx = abs(dx - dx.');
dy = abs(dy - dy.');
dz = max(dx, dy);

[Greg]

Ahh, I followed his link to a FEX submission and didn't dig any further.

[Walter Roberson]

I have seen several reports that implicit expansion is often slower than bsxfun. The details are likely to vary with release.

[Ansh]

Hi Greg, not sure if Jan has probably answered your question already but k is an input into the original function in the link in the original question. dx,dy and dz are already obtained from a different part of the code that Jan thankfully helped me to optimise with previously.

Answer 1

[Greg]

Taking some guesses:

eyeinds = eye(nObs,'logical');
dx2(eyeinds) = Inf;
dy2(eyeinds) = Inf;
dz2(eyeinds) = Inf;
% Edit to incorporate Jan's suggestions:
[EpsSample2, NnSample2] = mink(dz2,k,2);
Eps2 = EpsSample2(:,end);
Nn2 = NnSample2(:,end);
% [EpsSample2, NnSample2] = sort(dz2,2,'ascend');
% Eps2 = EpsSample2(:,k);
% Nn2 = NnSample2(:,k);
nx12 = sum(dx2 < Eps2,2);
ny12 = sum(dy2 < Eps2,2);
nx22 = sum(dx2 <= Eps2,2);
ny22 = sum(dy2 <= Eps2,2);

[Greg]

On my machine, with nObs = 10000;, I went from 21.5 to 12.5 seconds.

[Jan]

@Greg: Fine! +1

[Ansh]

As per my comment to Jan's answer below, I will download the newest version of MATLAB and test mink vs your unedited answer. Thanks for all the help from both of you, and hopefully I shall post soon with an update so I can accept an answer.

[Ansh]

Strange, I ran the code with the mink suggestion and X and Y = randn(4364,1) and k=3 (for definitions of X and Y see the link posted in the question and in Jan's link that he posted as a comment under the question). It should give a number quite close 0 (the mutual information between two independent time series), but with the mink code it doesn't give that at all? The code in Greg's answer but without mink and using sort gives the correct answer though?

[Jan]

What exactly does "it doesn't give that at all" mean? Please explain, what you do get. This would give us a chance to understand what happens. Post your code in addition.

These 2 codes should give the same output:

[EpsSample2, NnSample2] = sort(dz2, 2, 'ascend');
Eps2 = EpsSample2(:,k);
Nn2  = NnSample2(:,k);

[EpsSample2_2, NnSample2_2] = mink(dz2, 2, k);
Eps2_2 = EpsSample2_2(:, k);
Nn2_2  = NnSample2_2(:, k);

Do they? I cannot check this, because I do not have 2017b.

[Greg]

We both messed up. The arguments to mink are reversed.

mink(dz2,k,2);

Answer updated.

[Ansh]

Yes the updated answer is now working and is giving about 30% speed up from the previous modification. I'm very happy with the speed up of the code overall and am happy to accept this answer. Thanks Greg and Jan (again)!

[Greg]

Happy to help!

Answer 2

[Jan]

Sorting the complete array is a waste of time if all you need is the k.th smallest element. Under Matlab 2017b see mink to solve this much faster.

A further idea: Omit setting the diagonal to Inf, but leave it at zero. Then the result of nx12 etc. is 1 to large.

[Eps2, Nn2] = mink(dz2, 2, k);  
nx12 = sum(dx2 < Eps2, 2) - 1;
ny12 = sum(dy2 < Eps2, 2) - 1;
nx22 = sum(dx2 <= Eps2, 2) - 1;
ny22 = sum(dy2 <= Eps2, 2) - 1;

[Greg]

Leaving at zero only works if all data is guaranteed positive. And that throws off the min counting in dz. I didn't want to make that assumption.

Just noticed your version leaves Eps2 as a non-vector for k > 1. I've edited my answer to incorporate the awesome suggestion to use mink.

[Ansh]

Thanks for the input again Jan and to Greg too! Unfortunately I don't have MATLAB 2017b - I'm using MATLAB 2016b currently so I'll download the newest version (might as well I guess) and re test to see how much faster it is (hopefully the speed is up is significant).

And yes Greg luckily all the data is positive since dx, dy and dz are essentially distances.