Vectorization of nested for loops for following code

Can someone please help me to vectorize following code to avoid for loops. A similar kind of code is being called in one of my program too many times, hence need to optimize for efficiency purpose. Note that I don't want to use sum() or mean2() function as those are resulting efficiency reduction (at least when matrix size is 2x2). Thanks

Code :

    rangeDim = 4;
    domainDim = 2*rangeDim;
    dDomain = double(randi(10,domainDim));                        
    reducedDomain = double(zeros(rangeDim, rangeDim));
    for i=1:rangeDim
        i2 = i*2;
        for j=1:rangeDim
            j2 = j*2;
            reducedDomain(i,j) = 0.25*(dDomain(i2-1,j2-1)+dDomain(i2-1,j2)+dDomain(i2,j2-1)+dDomain(i2,j2));
            %reducedDomain(i,j) = mean2(dDomain(i2-1:i2,j2-1:j2));
            %reducedDomain(i,j) = 0.25*sum(sum(dDomain(i2-1:i2,j2-1:j2)));
        end
    end

[Matt J]

A similar kind of code is being called in one of my program too many times, hence need to optimize for efficiency purpose

Vectorizing operations on small chunks of data isn't going to be very beneficial. What you really want is to vectorize the outer loop, the one that performs this operation "many times".

Answer 1

[Guillaume]

By default, matlab creates double variables, therefore there's no need to convert zeros, randi, etc output to double. (and if it weren't it would be more efficient to specify the output type in the function call)

Your reduced domain is simply a part of the convolution of the original matrix with [1, 1;1, 1]. So:

rangeDim = 4;
domainDim = 2*rangeDim;
dDomain = randi(10,domainDim);  
reducedDomain = 0.25*conv2(dDomain, ones(2), 'valid');
reducedDomain = reducedDomain(1:2:end, 1:2:end)

[Matt J]

A shortcoming of this solution is that half of the computations done by the convolution get discarded in the final line,

reducedDomain = reducedDomain(1:2:end, 1:2:end)

[Guillaume]

Yes, another option would be to use blockproc from the image processing toolbox, but the cost of the anonymous function call may offset the gain in processing each block only once.

edit: Or as per the function link is your answer, you could split the matrix into cells and process each block that way. You've got the cost of a mat2cell, an anonymous function call and a for loop. You would have to benchmark to see if it is faster than computing a convolution and discarding half of it.

[Matt J]

Just to be clear, I'm advocating neither mat2cell, blockproc, nor conv2 for a situation like this. If you try sepblockfun, you should see that it gives an advantage over all of the above,

rangeDim = 2000;
domainDim = 2*rangeDim;
dDomain = randi(10,domainDim);  

tic;
 reducedDomain=sepblockfun(dDomain,[2,2],'mean');
toc
%Elapsed time is 0.071621 seconds.

tic;
reducedDomain = 0.25*conv2(dDomain, ones(2), 'valid');
reducedDomain = reducedDomain(1:2:end, 1:2:end);
toc
Elapsed time is 0.146615 seconds.%

tic
 e(1:rangeDim)=2;
 reducedDomain=cellfun(@(c) mean(c(:)), mat2cell(dDomain,e,e) );
toc
%Elapsed time is 76.570993 seconds.  

tic
 reducedDomain=blockproc(dDomain,[2,2],@(c) mean(c.data(:)) );
toc
%Elapsed time is 190.382548 seconds.

[Ashish KL]

Thanks for the help and all proposed solutions. I realized that sepblockfun() works with better efficiency for large size matrix, where as conv2() is performing better in average case.

Answer 2

[Matt J]

Using SEPBLOCKFUN from the File Exchange ( Download )

reducedDomain=sepblockfun(dDomain,[2,2],'mean');

or

reducedDomain=sepblockfun(dDomain,[2,2],'sum')/4;

Answer 3

[Alan Weiss]

1. Why all the double calls? It seems to me that they are completely superfluous.
2. MATLAB is a matrix language. I believe that you can find a matrix A so that reducedDomain = A*dDomain or something similar.
3. It is possible that conv2 could help.

Alan Weiss

MATLAB mathematical toolbox documentation

Answer 4

[Matt J]

I realized that sepblockfun() works with better efficiency for large size matrix, where as conv2() is performing better in average case.

That would be because sepblockfun has overhead from M-coded pre-checks and argument parsing that it performs in addition to the actual computation. However, you can customize the key lines from sepblockfun to this particular computation as below. I'm finding that this outperforms the conv2 approach even for rangeDim as small as 4.

    rangeDim = 4;
    domainDim = 2*rangeDim;
    dDomain = randi(10,domainDim);  

    N=10000;

    tic;
    for i=1:N
      X=reshape(dDomain,2,rangeDim,2,rangeDim);
      reducedDomain=reshape(sum(sum(X,1),3),rangeDim,rangeDim)/4;
    end
    toc
    %Elapsed time is 0.027305 seconds.   

    tic;
    for i=1:N
     reducedDomain = 0.25*conv2(dDomain, ones(2), 'valid');
     reducedDomain = reducedDomain(1:2:end, 1:2:end);
    end
    toc
    %Elapsed time is 0.050848 seconds.

[Ashish KL]

Hello Matt, Yes, this approach works better than conv2. Thanks for all the suggestions you have provided.