title: Efficient algorithm to calculate pairwise distances for a large set of locations in matlab
body: <p>I recently got a problem when calculating the spherical distance with large set of locations in matlab. Here is the problem:</p>
  
  <p>1.locations are latitudes and longitudes, stored in a vector locs with dimension nx2 where n = 10^5.</p>
  
  <p>2.I want to calculate the spherical distance between any two locations in the vector locs. Note that the resulted distance matrix may or may not sparse, since I cannot check it without distance. Let's assume that it's sparse, because we can truncate the distance with a number, say 1000. (The units are km), and the resulted distance matrix is enough for my usage.</p>
  
  <p>3.I tried two methods to calculate it, but each of them has drawbacks.</p>
  
  <p>PS: I run the code on Mac OS with MATLAB 2014 student version. The function great_circle_distance() is pretty like Matlab built-in function calculating the great circle distance between any two locations with geospatial coordinates. The problem here is not relevant to usage of this function.</p>
  
  <p>Thanks for all kinds of suggestions in advance. </p>
  
  <p>--Disc</p>
  
  <h2>Method1:</h2>
  
  <pre><code>dist = zeros(n, 1);
  dist_vec = [];
  dist_id = [];
  dist_id2 = [];
  for i =1:n
  dist = great_circle_distance(locs(i, :), locs);
  dist_temp=(dist&lt;300) &amp;(dist&gt;0);   % for example, consider the distance between 0 and 300
  [dist_id_temp, dist_id2_temp,dist_vtemp]=find(dist_temp);
  dist_vec_temp=dist(dist_id_temp);
  dist_id=[dist_id; dist_id_temp];
  dist_id2=[dist_id2; dist_id2_temp];
  dist_vec=[dist_vec; dist_vec_temp];
  
  if mod(i, 1000) == 0
  
      i
  
  end
  end
  dist_mat = sparse(dist_id, dist_id2, dist_vec);
  </code></pre>
  
  <p>The drawback here is that it took long time to finish this, at least 48 hours, but it consumed not too much memory around 2-5G.</p>
  
  <h2>Method2:</h2>
  
  <pre><code>d = pdist(locs, @great_circle_distance);
  [i, j] = find(tril(true(n), -1)); % extract index below main diagonal
  d = d';
  a = [i, j, d];
  con = d == 0;   % specify condition. 
  a(con, :) = [];   % delete rows satisfying condition con.
  clear i j d;
  i = a(:, 1);
  j = a(:, 2);
  d = a(:, 3);
  dist_mat = sparse(i, j, d);
  </code></pre>
  
  <p>The drawback of this method is that it consumed too much memory (exceeding 15G) when calculating d in the code.</p>