<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/performance-when-algorithmics-meets-mathematics/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://privefl.github.io/blog/performance-when-algorithmics-meets-mathematics/"> Florian Privé</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<section>
<p>In this post, I talk about performance through an efficient algorithm I developed for finding closest points on a map. This algorithm uses both concepts from mathematics and algorithmics.</p>
<div>
<h2>Problem to solve</h2>
<p>This problem comes from a <a rel="nofollow" href="https://stackoverflow.com/q/49863185/6103040" target="_blank">recent question on StackOverflow</a>.</p>
<blockquote>
<p>I have two matrices, one is 200K rows long, the other is 20K. For each row (which is a point) in the first matrix, I am trying to find which row (also a point) in the second matrix is closest to the point in the first matrix. This is the first method that I tried on a sample dataset:</p>
</blockquote>
<div>
<pre><code><span># Test dataset: longitude and latitude</span>
pixels.latlon &lt;-<span> </span><span>cbind</span>(<span>runif</span>(<span>200000</span>, <span>min =</span> <span>-</span><span>180</span>, <span>max =</span> <span>-</span><span>120</span>),
                       <span>runif</span>(<span>200000</span>, <span>min =</span> <span>50</span>, <span>max =</span> <span>85</span>))
grwl.latlon &lt;-<span> </span><span>cbind</span>(<span>runif</span>(<span>20000</span>, <span>min =</span> <span>-</span><span>180</span>, <span>max =</span> <span>-</span><span>120</span>),
                     <span>runif</span>(<span>20000</span>, <span>min =</span> <span>50</span>, <span>max =</span> <span>85</span>))
<span># Calculate the distance matrix</span>
<span>library</span>(geosphere)
dist.matrix &lt;-<span> </span><span>distm</span>(pixels.latlon, grwl.latlon, <span>fun =</span> distHaversine)
<span># Pick out the indices of the minimum distance</span>
rnum &lt;-<span> </span><span>apply</span>(dist.matrix, <span>1</span>, which.min)</code></pre>
</div>
<p>At first, this problem was a memory problem because <code>dist.matrix</code> would take 30GB.</p>
<p>A simple solution to overcome this memory problem has been proposed:</p>
<div>
<pre><code><span>library</span>(geosphere)
rnum &lt;-<span> </span><span>apply</span>(pixels.latlon, <span>1</span>, <span>function</span>(x) {
  dm &lt;-<span> </span><span>distm</span>(x, grwl.latlon, <span>fun =</span> distHaversine)
  <span>which.min</span>(dm)
})</code></pre>
</div>
<p>Yet, a second problem remains, this solution would take 30-40 min to run.</p>
</div>
<div>
<h2>First idea of improvement</h2>
<p>In the same spirit as with this <a rel="nofollow" href="https://adv-r.hadley.nz/profiling.html#t-test" target="_blank">case study in book <em>Advanced R</em></a>, let us see the source code of the <code>distHaversine</code> function and see if we can adapt it for our particular problem.</p>
<div>
<pre><code><span>library</span>(geosphere)
distHaversine</code></pre>
</div>
<pre><code>## function (p1, p2, r = 6378137) 
## {
##     toRad &lt;- pi/180
##     p1 &lt;- .pointsToMatrix(p1) * toRad
##     if (missing(p2)) {
##         p2 &lt;- p1[-1, ]
##         p1 &lt;- p1[-nrow(p1), ]
##     }
##     else {
##         p2 &lt;- .pointsToMatrix(p2) * toRad
##     }
##     p = cbind(p1[, 1], p1[, 2], p2[, 1], p2[, 2], as.vector(r))
##     dLat &lt;- p[, 4] - p[, 2]
##     dLon &lt;- p[, 3] - p[, 1]
##     a &lt;- sin(dLat/2) * sin(dLat/2) + cos(p[, 2]) * cos(p[, 4]) * 
##         sin(dLon/2) * sin(dLon/2)
##     a &lt;- pmin(a, 1)
##     dist &lt;- 2 * atan2(sqrt(a), sqrt(1 - a)) * p[, 5]
##     return(as.vector(dist))
## }
## </code></pre></div></section>