 <p>You can get rid of that kernel loop to calculate <code>D</code> with this <a href="http://in.mathworks.com/help/matlab/ref/bsxfun.html" rel="nofollow"><code>bsxfun</code></a> based <a href="http://in.mathworks.com/help/matlab/matlab_prog/vectorization.html" rel="nofollow"><code>vectorized</code></a> approach -</p> <pre><code>D = squeeze(sum(bsxfun(@min,A,permute(B,[3 2 1])),2)) </code></pre> <p>Or avoid <code>squeeze</code> with this modification -</p> <pre><code>D = sum(bsxfun(@min,permute(A,[1 3 2]),permute(B,[3 1 2])),3) </code></pre> <p>If the calculations of <code>D</code> involve <code>max</code> instead of <code>min</code>, just replace <code>@min</code> with <code>@max</code> there.</p> <hr> <p><strong>Explanation:</strong> The way <code>bsxfun</code> works is that it does <em>expansion</em> on singleton dimensions and performs the operation as listed with <code>@</code> inside its call. Now, this expansion is basically how one achieves vectorized solutions that replace for-loops. By <code>singleton dimensions</code> in arrays, we mean dimensions of <code>1</code> in them. </p> <p>In many cases, singleton dimensions aren't already present and for vectorization with <code>bsxfun</code>, we need to create <code>singleton dimensions</code>. One of the tools to do so is with <a href="http://in.mathworks.com/help/matlab/ref/permute.html" rel="nofollow"><code>permute</code></a>. That's basically all about the way vectorized approach stated earlier would work.</p> <p>Thus, your kernel code -</p> <pre><code>... case {'dist_histint','ih'} [m,d]=size(A); [n,d1]=size(B); if (d ~= d1) error('column length of A (%d) != column length of B (%d)\n',d,d1); end % With the MATLAB JIT compiler the trivial implementation turns out % to be the fastest, especially for large matrices. D = zeros(m,n); for i=1:m % m is number of samples of A if (0==mod(i,1000)) fprintf('.'); end for j=1:n % n is number of samples of B D(i,j) = sum(min([A(i,:);B(j,:)]));%./max(A(:,i),B(:,j))); end end </code></pre> <p>reduces to -</p> <pre><code>... case {'dist_histint','ih'} [m,d]=size(A); [n,d1]=size(B); if (d ~= d1) error('column length of A (%d) != column length of B (%d)\n',d,d1); end D = squeeze(sum(bsxfun(@min,A,permute(B,[3 2 1])),2)) %// OR D = sum(bsxfun(@min,permute(A,[1 3 2]),permute(B,[3 1 2])),3) </code></pre> <p>I am assuming the line: <code>if (0==mod(i,1000)) fprintf('.'); end</code> isn't important to the calculations as it does printing of some message. </p>