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

entropyfor loopmutual informationperformance

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!

Best Answer

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

Related Solutions

MATLAB: Speed up for loop in this code for calculating mutual information (maybe using GPU computing)

For a 1000x1000 matrix, this is 6 times faster already:

% Version 1:
n  = size(X, 1);
X = X.';
Y = Y.';
dx = zeros(n, n);
dy = zeros(n, n);
for j = 1:n
   Xj = X(:, j);
   Yj = Y(:, j);
   for i = j+1:n
      dx(i,j) = sqrt(sum(bsxfun(@minus, X(:, i), Xj) .^ 2));
      dy(i,j) = sqrt(sum(bsxfun(@minus, Y(:, i), Yj) .^ 2));
      dx(j,i) = dx(i,j);
      dy(j,i) = dy(i,j);
   end
end
dz = max(dx, dy);

The original function took 29.5 sec (R2016b, Core2Duo, Win7/64), and the cleaned version 5.2 sec.

Here the data are processed columnwise, which is much faster because neighboring elements are accessed much faster in the memory. Then the comparison my max() is done outside the loop. And finally the resulting matrix is symmetric and you can omit the computation of X(:,i) and X(:,j) if you have the results for X(:,j) and X(:,i) already.

I tried to vectorized the inner loop:

% Version 2:
n  = size(X, 1);
X = X.';
Y = Y.';
dx = zeros(n, n);
dy = zeros(n, n);
for j = 1:n
   dx(j+1:n,j) = sqrt(sum(bsxfun(@minus, X(:, j+1:n), X(:, j)) .^ 2, 1));
   dy(j+1:n,j) = sqrt(sum(bsxfun(@minus, Y(:, j+1:n), Y(:, j)) .^ 2, 1));
   dx(j,j+1:n) = dx(j+1:n,j);
   dy(j,j+1:n) = dy(j+1:n,j);
end
dz = max(dx, dy);

But this takes 21 sec for 1000x1000 arrays. But for smaller 100x100 inputs it is faster: 1.2 sec instead of 2.2 sec (100 iterations).

Now you have an efficient function to start a parallelization or computation on the GPU. Maybe this is useful (I cannot test it):

% Version 3:
parfor v = 1:2
  if v == 1
    for j = 1:n
      dx(j+1:n, j) = sqrt(sum((X(:, j+1:n) - X(:, j)) .^ 2, 1));
      dx(j, j+1:n) = dx(j+1:n, j);
    end
  else
    for j = 1:n
      dy(j+1:n, j) = sqrt(sum((Y(:, j+1:n) - Y(:, j)) .^ 2, 1));
      dy(j, j+1:n) = dy(j+1:n, j);
    end
  end
end

But parfor for the inner loop will use more cores.

MATLAB: How to reduce number of loops

Ok, well, you could adapt the following to your case then. I am not checking in-depth that it is working though, and there is room for improvement as I didn't optimize indexing for managing properly boundaries..

 n    = 100 ;
 sig  = 20 * rand(n,n,n) ; 
 LX   = size(sig,1) ;
 LY   = size(sig,2) ;
 LZ   = size(sig,3) ;
 dx   = 0.1 + rand(LX, 1) ; % Made irregular for testing..
 dy   = 0.8 + zeros(LY, 1) ;
 dz   = 0.3 + zeros(LZ, 1) ;
 Ctop = zeros(size(sig)) ;
 % - Your code.
 tic ;
 for i = 2:LX-1
    for j = 2:LY-1
        for k = 2:LZ-1
      Ctop(i,j,k) = (-1/dz(k-1))*((sig(i-1,j  ,k-1)*((dx(i-1)*dy(j  ))/4))+...
                                  (sig(i  ,j  ,k-1)*((dx(i  )*dy(j  ))/4))+...
                                  (sig(i-1,j-1,k-1)*((dx(i-1)*dy(j-1))/4))+...
                                  (sig(i  ,j-1,k-1)*((dx(i  )*dy(j-1))/4)));
       end
    end
 end
 toc
 % - Meshgrid based code.
 Ctop_mesh = zeros(size(Ctop)) ;
 tic ;
 [dy_m0, dx_m0] = meshgrid(dy(2:end-1), dx(2:end-1)) ;
 [dy_m1, dx_m1] = meshgrid(dy(1:end-2), dx(1:end-2)) ;
 sig_m00 = sig(2:end-1,2:end-1,:) ;
 sig_m10 = sig(1:end-2,2:end-1,:) ;
 sig_m01 = sig(2:end-1,1:end-2,:) ;
 sig_m11 = sig(1:end-2,1:end-2,:) ;
 for k = 2 : LZ-1
      Ctop_mesh(2:end-1,2:end-1,k) = (-1/dz(k-1)) / 4 * ...
          (sig_m10(:,:,k-1) .* dx_m1 .* dy_m0 + ...
           sig_m00(:,:,k-1) .* dx_m0 .* dy_m0 +...
           sig_m11(:,:,k-1) .* dx_m1 .* dy_m1 +...
           sig_m01(:,:,k-1) .* dx_m0 .* dy_m1) ;
 end
 toc
 % - Check.
 tmp = abs(Ctop_mesh(2:end-1,2:end-1,2:end-1)-Ctop(2:end-1,2:end-1,2:end-1)) ;
 all(tmp(:) < 1e-12)

This code outputs:

 Elapsed time is 6.281855 seconds.          % 3 nested loops
 Elapsed time is 0.408983 seconds.          % meshgrid
 ans =                                      % check
     1

Note that it is possible to get rid of the last FOR loop and work with a 3D array for dz built this way:

 dz_m1 = permute(repmat(dz(1:end-1),[1,LX-1,LY-1]), [2,3,1]) ;

but I am not sure that it is more efficient (especially if you want to PARFOR the remaining loop).

Best Answer

Related Solutions

MATLAB: Speed up for loop in this code for calculating mutual information (maybe using GPU computing)

MATLAB: How to reduce number of loops

Related Question