title: Efficiently calculating weighted distance in MATLAB
body: <p><a href="https://stackoverflow.com/questions/23911670/efficiently-compute-pairwise-squared-euclidean-distance-in-matlab">Several</a> <a href="https://stackoverflow.com/questions/25780633/calculating-euclidean-distance-of-pairs-of-3d-points-in-matlab">posts</a> <a href="https://statinfer.wordpress.com/2011/11/14/efficient-matlab-i-pairwise-distances/" rel="nofollow noreferrer">exist</a> about efficiently calculating pairwise distances in MATLAB. These posts tend to concern quickly calculating euclidean distance between large numbers of points.</p>
  
  <p>I need to create a function which quickly calculates the pairwise differences between smaller numbers of points (typically less than 1000 pairs). Within the grander scheme of the program i am writing, this function will be executed many thousands of times, so even small gains in efficiency are important. The function needs to be flexible in two ways:</p>
  
  <ol>
  <li>On any given call, the distance metric can be euclidean OR city-block.</li>
  <li>The dimensions of the data are weighted.</li>
  </ol>
  
  <p>As far as i can tell, no solution to this particular problem has been posted. The statstics toolbox offers <a href="http://www.mathworks.com/help/stats/pdist.html" rel="nofollow noreferrer">pdist</a> and <a href="http://www.mathworks.com/help/stats/pdist2.html" rel="nofollow noreferrer">pdist2</a>, which accept many different distance functions, but not weighting. I have seen extensions of these functions that allow for weighting, but these extensions do not allow users to select different distance functions. </p>
  
  <p>Ideally, i would like to avoid using functions from the statistics toolbox (i am not certain the user of the function will have access to those toolboxes).</p>
  
  <p>I have written two functions to accomplish this task. The first uses tricky calls to repmat and permute, and the second simply uses for-loops.</p>
  
  <pre><code>function [D] = pairdist1(A, B, wts, distancemetric)
  
  % get some information about the data
      numA = size(A,1);
      numB = size(B,1);
  
      if strcmp(distancemetric,'cityblock')
          r=1;
      elseif strcmp(distancemetric,'euclidean')
          r=2;
      else error('Function only accepts "cityblock" and "euclidean" distance')
      end
  
  %   format weights for multiplication
      wts = repmat(wts,[numA,1,numB]);
  
  %   get featural differences between A and B pairs
      A = repmat(A,[1 1 numB]);
      B = repmat(permute(B,[3,2,1]),[numA,1,1]);
      differences = abs(A-B).^r;
  
  %   weigh difference values before combining them
      differences = differences.*wts;
      differences = differences.^(1/r);
  
  %   combine features to get distance
      D = permute(sum(differences,2),[1,3,2]);
  end
  </code></pre>
  
  <p>AND:</p>
  
  <pre><code>function [D] = pairdist2(A, B, wts, distancemetric)
  
  % get some information about the data
      numA = size(A,1);
      numB = size(B,1);
  
      if strcmp(distancemetric,'cityblock')
          r=1;
      elseif strcmp(distancemetric,'euclidean')
          r=2;
      else error('Function only accepts "cityblock" and "euclidean" distance')
      end
  
  %   use for-loops to generate differences
      D = zeros(numA,numB);
      for i=1:numA
          for j=1:numB
              differences = abs(A(i,:) - B(j,:)).^(1/r);
              differences = differences.*wts;
              differences = differences.^(1/r);    
              D(i,j) = sum(differences,2);
          end
      end
  end
  </code></pre>
  
  <p>Here are the performance tests:</p>
  
  <pre><code>A = rand(10,3);
  B = rand(80,3);
  wts = [0.1 0.5 0.4];
  distancemetric = 'cityblock';
  
  
  tic
  D1 = pairdist1(A,B,wts,distancemetric);
  toc
  
  tic
  D2 = pairdist2(A,B,wts,distancemetric);
  toc
  
  Elapsed time is 0.000238 seconds.
  Elapsed time is 0.005350 seconds.
  </code></pre>
  
  <p>Its clear that the repmat-and-permute version works much more quickly than the double-for-loop version, at least for smaller datasets. But i also know that calls to repmat often slow things down, however. So I am wondering if anyone in the SO community has any advice to offer to improve the efficiency of either function!</p>
  
  <h2>EDIT</h2>
  
  <p>@Luis Mendo offered a nice cleanup of the repmat-and-permute function using <a href="http://www.mathworks.com/help/matlab/ref/bsxfun.html" rel="nofollow noreferrer">bsxfun</a>. I compared his function with my original on datasets of varying size:</p>
  
  <p><a href="https://i.stack.imgur.com/54pZS.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/54pZS.jpg" alt="comparison"></a></p>
  
  <p>As the data become larger, the bsxfun version becomes the clear winner!</p>
  
  <h2>EDIT #2</h2>
  
  <p>I have finished writing the function and it is available on github [<a href="https://github.com/noconaway/pairdist" rel="nofollow noreferrer">link</a>]. I ended up finding a pretty good vectorized method for computing euclidean distance [<a href="http://trunghuyduong.blogspot.com/2011/09/pair-wise-weighted-euclidean-distance.html" rel="nofollow noreferrer">link</a>], so i use that method in the euclidean case, and i took @Divakar's <a href="https://stackoverflow.com/a/33247023/3521179">advice</a> for city-block. It is still not as fast as pdist2, but its must faster than either of the approaches i laid out earlier in this post, and easily accepts weightings.</p>