how to make the .m code faster ?

I wrote a matlab code but the execution time is more than 25 minutes . I think it is because the for loop and the index. so i have a matrix of data H(435*4258) , and i am trying to take vaulues from it depend on another values from another matrix. x,y have the same size (435*4258) like H and contain the location of the data in H . the output matrix has c(869*869) . what Iam trying to do is making a rectangle around every point and search for each point in the matrix data is in this range and take the average . can somebody help me to make it faster and less than 2 seconds. the code is like interpolation. is there any way to make it faster ? the code :

 destance = 7.5;
d = distance/2;
m = 869;
n = 869;
rx = 3255;
c= zeros(m ,n);
newx = -rx:distance:+rx;
newy = -rx:distance:+rx;
for i = 1:m
    yc = newy(i)+d;
    yf = newy(i)-d;
    for j = 1:n
          xc = newx(j)+d;
          xf = newx(j)-d;
         ind = find(x>= xf & x<=xc & y >=yf & y<=yc); 
         if isempty(ind)
            c(i,j) = NaN; 
           else
           p = mean(h(ind));
           c(i,j) = p;
         end
    end
end

[Adam]

Look at and use

doc profile

before you start making any assumptions of exactly which part of code is slow. It is a very simple and easy to use profiler compared to many you get for other languages.

[Guillaume]

The code is bound to be slow. For each i, it rescans the x array which does not change between i for a given j. For each j, it rescans the y array which does not change for a given i. In other words, it performs m*n scans where only m+n scans are required at most.

[Stephen23]

https://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html

https://www.mathworks.com/help/matlab/matlab_prog/preallocating-arrays.html

https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html

https://www.mathworks.com/matlabcentral/answers/228557-experts-of-matlab-how-did-you-learn-any-advice-for-beginner-intermediate-users

Answer 1

[Guillaume]

First, a piece of advice, don't hardcode values that matlab can easily calculate. In your code m and n are the number of elements in vectors newx and newy, so you shouldn't hardcode them. One day, you'll decide to change distance or rx and your program will fail because you'll forget to recalculate m and n. Much safer:

m = numel(newx); 
n = numel(newy);

Anyway, if I understood correctly, this will produce the same result as your code:

distance = 7.5;
rx = 3255;
newx = -rx:distance:rx;
newy = -rx:distance:ry;    %shouldn't there be a ry?

binx = discretize(x(:), newx - distance/2);  %find index of x bin centered on newx with width distance
biny = discretize(y(:), newy - distance/2);  %find index of y bin centered on newy with width distance
c = accumarray([biny, binx], h(:), [], @mean);  %calculate the mean of all the values that fall within a bin

And should be much faster than your code which scan over and over the x and y arrays.

[Jony Muller]

thanks for your answer but i got , Undefined function or variable 'accummaray'. also incase i want to replace nan members in c with the nearst value in H what shoud I do. mabey you mean accumarray but i got this Error using accumarray First input SUBS must contain positive integer subscripts. regards

[Guillaume]

Yes, typo, it should have been accumarray.

The error you get about SUBS is most likely because you have some coordinates outside your grid (resulting in NaNs in the output of discretize). Not sure what you want to do about that, to ignore these points:

binx = discretize(x(:), newx - distance/2);  %find index of x bin centered on newx with width distance
biny = discretize(y(:), newy - distance/2);  %find index of y bin centered on newy with width distance
ingrid = ~isnan(binx) & ~isnan(biny);
c = accumarray([biny(ingrid), binx(ingrid)], h(ingrid), [], @mean);  %calculate the mean of all the values that fall within a bin

[Jony Muller]

thanks super it works in 3 seconds (Elapsed time is 2.968575 seconds.) but if i do not want 0 in the out put ? mabey because i have zeros in mu input . you have an idea ? also if it is possible to make all less than 2 seconds. thanks

[Guillaume]

The fifth value of accumarray specifies the fill value for the grid square that are empty, so:

c = accumarray([biny(ingrid), binx(ingrid)], h(ingrid), [], @mean , NaN);

Short of writing mex code, I doubt you can make this any faster.

[Jony Muller]

thanks.

[Guillaume]

I made a silly mistake in the code I wrote in the comment to Jan's answer. Initialising the sum matrix to NaN is never going to sum to anything but NaN. The correct code is:

edited: fixed off by 1 calculation of size of bin matrices.

outsize = [2*rx/distance, 2*rx/distance];
binsum = zeros(outsize);
bincount = zeros(outsize);
for idx = 1:numel(h)
   binx = 1 + floor((x(idx)+rx)/distance);  %discretisation of x
   biny = 1 + floor((y(idx)+rx)/distance);  %discretisation of y   
   if binx > 0 & biny > 0 & binx <= outsize(1) & biny <= outsize(2)
      binsum(biny, binx)  = binsum(biny, binx) + h(idx);
      bincount(biny, binx) = bincount(biny, binx) + 1;
   end
end
binmean = binsum ./ bincount; %internally loops over the bins.

To my surprise it is much much faster than my original solution.

Be aware that there is a minor difference between the two solutions, the grid in this latest solution is centered on the mid-points between the -rx:distance:rx whereas your original code and my accumarray code center the grid on the -rx:distance:rx points.

I'm going to assume it was a mistake in your original code (and this was also why you got NaNs in the discretize calls)

[Jony Muller]

thanks it workes faster but please can you explain it to me step by step and why you did that because i did not understand it . thanks

[Guillaume]

"i did not understand it" Really? You can't have spent much time trying to understand it, particularly considering the lengthy explanation I gave in Jan's answer thread.

The code steps simultaneous through the <x, y, h> triplets. At each step, finds out which bin the x and y fall in, add h to that bin sum and add 1 to the bin count. When the loop is finished, divide the sum by the count to get the mean. As a bonus, 0/0 is NaN, so you automatically get NaN when nothing goes into the bin.

[Jony Muller]

thank you very much but the last thing i need is in the second selution i have NaN so if I want to replace this NaN with the nearst Value in H what should I do ? bincount (biny,binx) = min(sqrt((x-newx(binx)).*(x-newx(binx))+(y-newx(biny)).*(y-newx(biny))));

[Guillaume]

We spent all this time trying to get this to work as fast as possible and you want to kill all that performance by searching for the nearest value. That does not make much sense to me.

That line of code you wrote also makes no sense. I don't even understand why you'd try to put the minimum of something into the bin count.

This may work (untested):

[binycentres, binxcentres] = ndgrid((-rx:distance:rx) + distance/2);
nanbins = isnan(binmean);
disttocentre = hypot(binxcentres(nanbins)' - x(:), binycentres(nanbins)' - y(:));  %requires R2016b or later
[~, nearestidx] = min(disttocentre);
binmean(nanbins) = h(nearestidx);

[Jony Muller]

the problem i don't want to lose data and if i have NaN in the output I have to take the nearst value in H. I tried this but it doesn't work mabey because I have R2015b the error is "Error using - Matrix dimensions must agree." . I tried to but minimun because i have alot of zeros in bin count abd we divide by bin count also the number of zeros is thae same number og NaN

[Guillaume]

In version < R2016b:

disttocentre = hypot(bsxfun(@minus, binxcentres(nanbins)', x(:)), bsxfun(@minus, binycentres(nanbins)', y(:)));

Finding the nearest point for each nan bin is never going to be fast. You have to calculate the distance of each bin centre to each of your 1.8 million points. So that's 1.8 million * number of nan bins hypot calculations.

[Jony Muller]

Error using bsxfun Requested 1852230x163817 (2260.7GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or preference panel for more information.

[Guillaume]

Well use a loop then. More than 20% of your bins are NaN. That's a lot. as a result you're trying to do 303,426,761,910 hypot calculation at once, which is going to need a lot of memory.

[binycentres, binxcentres] = ndgrid((-rx:distance:rx) + distance/2);
for binidx = find(isnan(binmean))'
   binidx
   disttocentre = hypot(x(:) - binxcentres(binidx), y(:) - binycentres(binidx));
   [~, nearestidx] = min(disttocentre);
   binmean(binidx) = h(nearestidx);
end

That loop only does 1,852,230 hypot at once, but does it 163,817 times. You're trading speed for memory.

No matter what it's going to be very slow. You probably want to rethink what you're doing if you want to do this in real time.

Answer 2

[Jan]

I assume Guillaume's suggestion is faster. For a comparison thry this cleaned loop:

distance = 7.5;  % Not "destance"
d  = distance/2;
m  = 869;
n  = 869;
rx = 3255;
c  = nan(m, n);
newx = -rx:distance:+rx;
newy = -rx:distance:+rx;
for i = 1:m
  yc = newy(i) + d;
  yf = newy(i) - d;
  yi = (y >= yf & y <= yc);  % Once per loop only
  for j = 1:n
     xc  = newx(j) + d;
     xf  = newx(j) - d;
     ind = (x >= xf & x <= xc & yi); 
     if any(ind)
       c(i, j) = sum(h(ind)) / sum(ind);
     end
  end
end

@Jony Muller: Please run a TIC/TOC and post the results for the two methods. Thanks.

[Guillaume]

Really, the two loops should be separate since x and y are completely independent. You've taken out the multi-scanning of the y array but not the x.

[Jony Muller]

thanks for your answer also i got Error using / Matrix dimensions must agree. and if i want to replace nan in c with the nearst value in h what shoud I do ?

[Jony Muller]

Elapsed time is 1357.365452 seconds. and i change the sum(h(ind)) / sum(ind); to mean but i got every thing nan so no data > any idea ?

[Jan]

@Guillaume: I don't get your point.

@Jony: mean() has more overhead than sum. If you provide some meaningful input data, we could do some time measurements by out own. Optimizations are much easier, when e.g. the output of the profiler can be used. This code needs 1350 seconds?! This sounds really strange.

Nevertheless, if Guillaume's version works, you should use it in every case and comparing it with leaner loops is of accademic interest only.

[Jony Muller]

Guillaume's version works very fine but the time is 3 seconds what I need is less than 2 seconds because I want to make video and have every 2 seconds I have an image also I change sum to mean because I got an error with / after few minutes of excution ( Error using / Matrix dimensions must agree). and the output is nan so I have nothing.

[Jan]

Without having the corresponding inputs, I cannot guess, what happens. How large are the inputs? Even 3 seconds seems to be surprisingly slow.

I do not think, that the shown code can create a "/ Matrix dimensions must agree" error. Did you modify the code?

It would be very useful to have some data to play with. Then answering is not based on guessing.

[Jony Muller]

the Input is matrix of data H(435*4258) this data and conatins positive data also it could contain zeros.newx equal to newy with width (435*2 -1). the problem with "/" is when the function find the index and try to make the average i got this error. I did just main intead of sum/sum because i got an error and when the function finsh the time was 1357.365452 seconds with output without data just NaN. I tried to put sum./sum so I got (Assignment has more non-singleton rhs dimensions than non-singleton subscripts) I think the problem is with the two for loop it make every time check all elements of x and y and h. do you have any idea ? thx

[Jan]

@Jony: Again, please provide some input data. If you explain how the input data might look like, I would have to sit down and invent some Matlab code to produce some arrays which might match yours. But this would be a waste of my time, and even your time, when my speculations do not match your data.

I cannot imagine or guess, why mean could reply anything different from sum(h(ind)) / sum(ind). As long as ind is a vector, both epressions should be scalars.

I really want to help you. But without meaningful data I do not see a way do to this efficiently. So please post either the original data, or if it is enough to reproduce the problem some rand calls to create pseudo data. The sizes and values of x and y matters, perhaps they are sorted already. H is not huge, therefore I do not have the faintest idea what the computer is doing in the 1357 seconds.

[Jony Muller]

you are right but i have the input data for an image how do sugeest to send to show this data let us say a random positive data with zeros (not sorted) and also i am surprised why i have this error and with mean i got nothing with my original code i have outdata with nan and values but the time is 1860 seconds.

[Jan]

Can you attach the data as MAT files?

[Jony Muller]

Cannot attach this file because:

File size exceeds 5 MB. Try compressing the image file in an archive format like zip. the size is 120 MB .

[Jony Muller]

[Guillaume]

@Jan, "Guillaume: I don't get your point."

You've extracted yi = (y >= yf & y <= yc) out of the inner loop because it kept being recalculated for each x position. However, the same happens with x. The result x >= xf & x <= xc is the same for a given j regardless of the value of i, yet you recalculate it at each step of the i loop.

That is my point, the binning of the x is independent of the binning of the y. Therefore, it would make sense to have the two loops independent. Practically, with your approach, that's not possible (without a third loop).

[Guillaume]

@Jony,

Internally, my solution involves 4 loops (I assume, if discretize and accumarray are implemented optimally)

• one loop over the x values, for the first discretize
• one loop over the y values, for the second discretize
• one loop over corresponding x,y and h values by accumarray to accumulate the h values into their corresponding bin
• one loop over the bins by accumarray to calculate the bin mean

All these loops are implemented in mex files or directly within the matlab compiled code so they'll run much faster than anything you'd write in m code.

You could reduce it to two loops by doing the discretisation and accumulation in the same loop. Code in m would be something like:

binsum = nan(2*rx/distance + 1);
bincount = zeros(size(binsum));
for idx = 1:numel(h)
   binx = 1 + floor((x(idx)+rx)/distance);  %discretisation of x
   biny = 1 + floor((y(idx)+rx)/distance);  %discretisation of y   
   if binx > 0 & biny > 0 & binx <= m & biny <= m
      binsum(biny, binx)  = binsum(biny, binx) + h(idx);
      bincount(biny, binx) = bincount(biny, binx) + 1;
   end
end
binmean = binsum ./ bincount; %internally loops over the bins.

This is optimal in the number of loops. You could make it even faster, if you got rid of the branching if by ensuring that all x and y are within -rx:+rx.

However, unless implemented in mex, I doubt it'll be faster than my discretize + accumarray solution.

[Jan]

@Guillaume: Thanks, it is clear now. I only simplified the original loops in a trivial way. I'm afraid that x and y are such huge that sorting the complete matrix might exhaust the memory.

[Jony Muller]

@Guillaume @Jan Simon so what is the final solution. because in this I got NaN

binsum = nan(2*rx/distance + 1);
bincount = zeros(size(binsum));
for idx = 1:numel(h)
   binx = 1 + floor((x(idx)+rx)/distance);  %discretisation of x
   biny = 1 + floor((y(idx)+rx)/distance);  %discretisation of y   
   if binx > 0 & biny > 0 & binx <= m & biny <= m
      binsum(biny, binx)  = binsum(biny, binx) + h(idx);
      bincount(biny, binx) = bincount(biny, binx) + 1;
   end
end
binmean = binsum ./ bincount; %internally loops over the bins.

thanks for your help

[Guillaume]

"because in this I got NaN"

Yes, but you got it very fast ;) See as a comment to my answer for the fix.

[Jan]

@Jony: As long as you cannot provide input data, I cannot test some ideas I have. I understand that it is not trivial due to the file size. But this is your problem, so it is your turn to find a way to allow us to test suggestions and improvements.

But as far as I understand, Guillaume's code solves the problem already. Or does it reply something different from your code?

[Guillaume]

@Jan,

I tested with this:

x = rx*(2*rand(435, 4258) - 1);
y = rx*(2*rand(435, 4258) - 1);
h = randi(200, 435, 4258);

which corresponds more or less to tony's description.

I don't know why the data is arrange in matrix form. As far as I understand this has no significance. The three variables could be just vectors.