title: Compute the Cartesian Distance between points in high-dimensional space optimally
body: <p><strong>Summary:</strong></p>
  
  <p>My class SquareDistance computes the square of the Cartesian distance in four ways using methods with these names:</p>
  
  <ol>
  <li>Signed</li>
  <li>UnsignedBranching</li>
  <li>UnsignedDistribute</li>
  <li>CastToSignedLong</li>
  </ol>
  
  <p>The first one is fastest and uses signed integers, but my data must be unsigned (for reasons given below). The other three methods start with unsigned numbers. My goal is to write a method like those in SquareDistance that takes unsigned data and performs better than the three I already wrote, as close as possible in performance to #1. Code with benchmark results follows.  (unsafe code is permitted, if you think it will help.)</p>
  
  <p><strong>Details:</strong></p>
  
  <p>I am developing an algorithm to solve K-nearest neighbor problems using an index derived from the Hilbert curve. The execution time for the naive, linear scan algorithm grows in time quadratically with the number of points and linearly with the number of dimensions, and it spends all its time computing and comparing Cartesian distances. </p>
  
  <p>The motivation behind the special Hilbert index is to reduce the number of times that the distance function is called. However, it must still be called millions of times, so I must make it as fast as possible. (It is the most frequently called function in the program. A recent failed attempt to optimize the distance function doubled the program execution time from seven minutes to fifteen minutes, so no, this is not a premature or superfluous optimization.)</p>
  
  <p><strong>Dimensions</strong>: The points may have anywhere from ten to five thousand dimensions.</p>
  
  <p><strong>Constraints</strong>. I have two annoying constraints:</p>
  
  <ol>
  <li><p>The Hilbert transformation logic requires that the points be expressed as uint (unsigned integer) arrays. (The code was written by another, is magic and uses shifts, ANDs, ORs, and the like and can't be changed.) Storing my points as signed integers and incessantly casting them to uint arrays produced wretched performance, so I must at the very least store a uint array copy of each point.</p></li>
  <li><p>To improve efficiency, I made a signed integer copy of each point to speed up the distance calculations. This worked very well, but once I get to about 3,000 dimensions, I run out of memory! </p></li>
  </ol>
  
  <p>To save on memory, I removed the memoized signed integer arrays and tried to write an unsigned version of the distance calculation. My best results are 2.25 times worse than the signed integer version.</p>
  
  <p>The benchmarks create 1000 random points of 1000 dimensions each and perform distance calculations between every point and every other point, for 1,000,000 comparisons. Since I only care about the relative distance, I save time by not performing the square root.</p>
  
  <p>In debug mode:</p>
  
  <blockquote>
  <pre><code>SignedBenchmark                  Ratio: 1.000 Seconds: 3.739
  UnsignedBranchingBenchmark       Ratio: 2.731 Seconds: 10.212
  UnsignedDistributeBenchmark      Ratio: 3.294 Seconds: 12.320
  CastToSignedLongBenchmark        Ratio: 3.265 Seconds: 12.211
  </code></pre>
  </blockquote>
  
  <p>In release mode:</p>
  
  <pre><code> SignedBenchmark                  Ratio: 1.000 Seconds: 3.494
   UnsignedBranchingBenchmark       Ratio: 2.672 Seconds: 9.334
   UnsignedDistributeBenchmark      Ratio: 3.336 Seconds: 11.657
   CastToSignedLongBenchmark        Ratio: 3.471 Seconds: 12.127
  </code></pre>
  
  <p>The above benchmarks were run on a Dell with an Intel Core i7-4800MQ CPU @ 2.70GHz with 16 GB memory. My larger algorithm already uses the Task Parallel library for larger tasks, so it is fruitless to parallelize this inner loop.</p>
  
  <p><strong>Question:</strong> Can anyone think of a faster algorithm than UnsignedBranching?</p>
  
  <p>Below is my benchmark code.</p>
  
  <p><strong>UPDATE</strong></p>
  
  <p>This uses loop unrolling (thanks to @dasblinkenlight) and is 2.7 times faster:</p>
  
  <pre><code>public static long UnsignedLoopUnrolledBranching(uint[] x, uint[] y)
  {
      var distance = 0UL;
      var leftovers = x.Length % 4;
      var dimensions = x.Length;
      var roundDimensions = dimensions - leftovers;
  
      for (var i = 0; i &lt; roundDimensions; i += 4)
      {
          var x1 = x[i];
          var y1 = y[i];
          var x2 = x[i+1];
          var y2 = y[i+1];
          var x3 = x[i+2];
          var y3 = y[i+2];
          var x4 = x[i+3];
          var y4 = y[i+3];
          var delta1 = x1 &gt; y1 ? x1 - y1 : y1 - x1;
          var delta2 = x2 &gt; y2 ? x2 - y2 : y2 - x2;
          var delta3 = x3 &gt; y3 ? x3 - y3 : y3 - x3;
          var delta4 = x4 &gt; y4 ? x4 - y4 : y4 - x4;
          distance += delta1 * delta1 + delta2 * delta2 + delta3 * delta3 + delta4 * delta4;
      }
      for (var i = roundDimensions; i &lt; dimensions; i++)
      {
          var xi = x[i];
          var yi = y[i];
          var delta = xi &gt; yi ? xi - yi : yi - xi;
          distance += delta * delta;
      }
      return (long)distance;
  }
  </code></pre>
  
  <p><strong>SquareDistance.cs:</strong></p>
  
  <pre><code>using System;
  using System.Collections.Generic;
  using System.Linq;
  using System.Text;
  using System.Threading.Tasks;
  
  namespace DistanceBenchmark
  {
      /// &lt;summary&gt;
      /// Provide several alternate methods for computing the square of the Cartesian distance
      /// to allow study of their relative performance.
      /// &lt;/summary&gt;
      public static class SquareDistance
      {
          /// &lt;summary&gt;
          /// Compute the square of the Cartesian distance between two N-dimensional points
          /// with calculations done on signed numbers using signed arithmetic, 
          /// a single multiplication and no branching.
          /// &lt;/summary&gt;
          /// &lt;param name="x"&gt;First point.&lt;/param&gt;
          /// &lt;param name="y"&gt;Second point.&lt;/param&gt;
          /// &lt;returns&gt;Square of the distance.&lt;/returns&gt;
          public static long Signed(int[] x, int[] y)
          {
              var distance = 0L;
              var dimensions = x.Length;
              for (var i = 0; i &lt; dimensions; i++)
              {
                  var delta = x[i] - y[i];
                  distance += delta * delta;
              }
              return distance;
          }
  
          /// &lt;summary&gt;
          /// Compute the square of the Cartesian distance between two N-dimensional points
          /// with calculations done on unsigned numbers using unsigned arithmetic, a single multiplication
          /// and a branching instruction (the ternary operator).
          /// &lt;/summary&gt;
          /// &lt;param name="x"&gt;First point.&lt;/param&gt;
          /// &lt;param name="y"&gt;Second point.&lt;/param&gt;
          /// &lt;returns&gt;Square of the distance.&lt;/returns&gt;
          public static long UnsignedBranching(uint[] x, uint[] y)
          {
              var distance = 0UL;
              var dimensions = x.Length;
              for (var i = 0; i &lt; dimensions; i++)
              {
                  var xi = x[i];
                  var yi = y[i];
                  var delta = xi &gt; yi ? xi - yi : yi - xi;
                  distance += delta * delta;
              }
              return (long)distance;
          }
  
          /// &lt;summary&gt;
          /// Compute the square of the Cartesian distance between two N-dimensional points
          /// with calculations done on unsigned numbers using unsigned arithmetic and the distributive law,
          /// which requires four multiplications and no branching.
          /// 
          /// To prevent overflow, the coordinates are cast to ulongs.
          /// &lt;/summary&gt;
          /// &lt;param name="x"&gt;First point.&lt;/param&gt;
          /// &lt;param name="y"&gt;Second point.&lt;/param&gt;
          /// &lt;returns&gt;Square of the distance.&lt;/returns&gt;
          public static long UnsignedDistribute(uint[] x, uint[] y)
          {
              var distance = 0UL;
              var dimensions = x.Length;
              for (var i = 0; i &lt; dimensions; i++)
              {
                  ulong xi = x[i];
                  ulong yi = y[i];
                  distance += xi * xi + yi * yi - 2 * xi * yi;
              }
              return (long)distance;
          }
  
          /// &lt;summary&gt;
          /// Compute the square of the Cartesian distance between two N-dimensional points
          /// with calculations done on unsigned numbers using signed arithmetic, 
          /// by first casting the values into longs.
          /// &lt;/summary&gt;
          /// &lt;param name="x"&gt;First point.&lt;/param&gt;
          /// &lt;param name="y"&gt;Second point.&lt;/param&gt;
          /// &lt;returns&gt;Square of the distance.&lt;/returns&gt;
          public static long CastToSignedLong(uint[] x, uint[] y)
          {
              var distance = 0L;
              var dimensions = x.Length;
              for (var i = 0; i &lt; dimensions; i++)
              {
                  var delta = (long)x[i] - (long)y[i];
                  distance += delta * delta;
              }
              return distance;
          }
  
      }
  }
  </code></pre>
  
  <p><strong>RandomPointFactory.cs:</strong></p>
  
  <pre><code>using System;
  using System.Collections.Generic;
  using System.Linq;
  using System.Text;
  using System.Threading.Tasks;
  
  namespace DistanceBenchmark
  {
      public static class RandomPointFactory
      {
          /// &lt;summary&gt;
          /// Get a random list of signed integer points with the given number of dimensions to use as test data.
          /// &lt;/summary&gt;
          /// &lt;param name="recordCount"&gt;Number of points to get.&lt;/param&gt;
          /// &lt;param name="dimensions"&gt;Number of dimensions per point.&lt;/param&gt;
          /// &lt;returns&gt;Signed integer test data.&lt;/returns&gt;
          public static IList&lt;int[]&gt; GetSignedTestPoints(int recordCount, int dimensions)
          {
              var testData = new List&lt;int[]&gt;();
              var random = new Random(DateTime.Now.Millisecond);
  
              for (var iRecord = 0; iRecord &lt; recordCount; iRecord++)
              {
                  int[] point;
                  testData.Add(point = new int[dimensions]);
                  for (var d = 0; d &lt; dimensions; d++)
                      point[d] = random.Next(100000);
              }
              return testData;
          }
  
          /// &lt;summary&gt;
          /// Get a random list of unsigned integer points with the given number of dimensions to use as test data.
          /// &lt;/summary&gt;
          /// &lt;param name="recordCount"&gt;Number of points to get.&lt;/param&gt;
          /// &lt;param name="dimensions"&gt;Number of dimensions per point.&lt;/param&gt;
          /// &lt;returns&gt;Unsigned integer test data.&lt;/returns&gt;
          public static IList&lt;uint[]&gt; GetUnsignedTestPoints(int recordCount, int dimensions)
          {
              var testData = new List&lt;uint[]&gt;();
              var random = new Random(DateTime.Now.Millisecond);
  
              for (var iRecord = 0; iRecord &lt; recordCount; iRecord++)
              {
                  uint[] point;
                  testData.Add(point = new uint[dimensions]);
                  for (var d = 0; d &lt; dimensions; d++)
                      point[d] = (uint)random.Next(100000);
              }
              return testData;
          }
      }
  }
  </code></pre>
  
  <p><strong>Program.cs:</strong></p>
  
  <pre><code>using System;
  using System.Collections.Generic;
  using System.Diagnostics;
  using System.Linq;
  using System.Text;
  using System.Threading.Tasks;
  
  namespace DistanceBenchmark
  {
      public class Program
      {
          private static IList&lt;int[]&gt; SignedTestData = RandomPointFactory.GetSignedTestPoints(1000, 1000);
          private static IList&lt;uint[]&gt; UnsignedTestData = RandomPointFactory.GetUnsignedTestPoints(1000, 1000);
  
          static void Main(string[] args)
          {
              var baseline = TimeIt("SignedBenchmark", SignedBenchmark);
              TimeIt("UnsignedBranchingBenchmark", UnsignedBranchingBenchmark, baseline);
              TimeIt("UnsignedDistributeBenchmark", UnsignedDistributeBenchmark, baseline);
              TimeIt("CastToSignedLongBenchmark", CastToSignedLongBenchmark, baseline);
              TimeIt("SignedBenchmark", SignedBenchmark, baseline);
              Console.WriteLine("Done. Type any key to exit.");
              Console.ReadLine();
          }
  
          public static void SignedBenchmark()
          {
              foreach(var p1 in SignedTestData)
                  foreach (var p2 in SignedTestData)
                      SquareDistance.Signed(p1, p2);
          }
  
          public static void UnsignedBranchingBenchmark()
          {
              foreach (var p1 in UnsignedTestData)
                  foreach (var p2 in UnsignedTestData)
                      SquareDistance.UnsignedBranching(p1, p2);
          }
  
          public static void UnsignedDistributeBenchmark()
          {
              foreach (var p1 in UnsignedTestData)
                  foreach (var p2 in UnsignedTestData)
                      SquareDistance.UnsignedDistribute(p1, p2);
          }
  
          public static void CastToSignedLongBenchmark()
          {
              foreach (var p1 in UnsignedTestData)
                  foreach (var p2 in UnsignedTestData)
                      SquareDistance.CastToSignedLong(p1, p2);
          }
  
          public static double TimeIt(String testName, Action benchmark, double baseline = 0.0)
          {
              var stopwatch = new Stopwatch();
              stopwatch.Start();
              benchmark();
              stopwatch.Stop();
              var seconds = stopwatch.Elapsed.TotalSeconds;
              var ratio = baseline &lt;= 0 ? 1.0 : seconds/baseline;
              Console.WriteLine(String.Format("{0,-32} Ratio: {1:0.000} Seconds: {2:0.000}", testName, ratio, seconds));
              return seconds;
          }
      }
  }
  </code></pre>