Repository: scipy/scipy
Pull request: ENH: Add fast path for minkowski distance with p=1,2 and support p=inf
I noticed that `spatial.distance_matrix(x, y)` was much faster than `spatial.distance.cdist(x, y, metric='minkowski')`. But only when `p` is left at the default of two. I'm guessing NumPy special cases `x**2` and turns it into `x*x`.

This special cases `p=2` and `p=1` in `cdist` and `pdist` so that it doesn't use `pow`. With this, `cdist` is 10x faster than `distance_matrix` for `p=1`, `p=2` and also marginally faster for other values of `p` as well.

Edit: I've included a fix for `p=inf` case as well (fixes gh-8776) since the fix is in a very similar vein.

Problem statement:
cdist behavior with Minkowsky and np.inf
Using `scipy.spatial.distance.cdist` with the Minkowski distance and `p=np.inf` does not produce expected result. In particular, I'm not sure how `cdist` is producing all distance of only 1's (see Minkowski below). Using Chebyshev distance produce expected results, but they should be equivalent for $p=\infty$.

```
import numpy as np
import scipy.spatial as sp

a = np.array([[1.11132, 1.15693], [0.98214, 1.01018]])
b = np.array([[1.14994, 1.15386 ], [0.973258, 1.17015 ], [0.939203, 1.09205 ]])

print("a:\n{}".format(a))
print("b:\n{}".format(b))

mink = sp.distance.cdist(a, b, 'minkowski', p=np.inf)
print("Minkowski:\n{}".format(mink))

cheb = sp.distance.cdist(a, b, 'chebyshev')
print("Chebyshev:\n{}".format(cheb))

a:
[[1.11132 1.15693]
 [0.98214 1.01018]]
b:
[[1.14994  1.15386 ]
 [0.973258 1.17015 ]
 [0.939203 1.09205 ]]
Minkowski:
[[1. 1. 1.]
 [1. 1. 1.]]
Chebyshev:
[[0.03862  0.138062 0.172117]
 [0.1678   0.15997  0.08187 ]]
```

### Scipy/Numpy/Python version information:
```
scipy.__version__  = 1.0.0
numpy.__version__ = 1.14.2
```

Thanks,
Jesse

Hints from issue discussion:
The reason it returns incorrect results are because the implemntation doesn't special-case ``p=np.inf``, and so ends up raising every distance to the zero power; see the implementation here:

https://github.com/scipy/scipy/blob/5e22b2e447cec5588fb42303a1ae796ab2bf852d/scipy/spatial/src/distance_impl.h#L338-L349

Two possible fixes:

1. Change ``minkowski_distance`` so that it falls back to ``chebyshev_distance_double`` when ``p`` is infinite

2. Raise a ValueError when p is infinite and tell the user to use chebyshev distance instead
I'm happy to work on either fix... any thoughts on the preferred approach?
@jakevdp seems that the first approach would be the least disruptive. I can't think of any reason to raise an error since falling back to `chebyshev_distance_double` provides the user with the desired answer for infinite `p`.

Commit message: ENH: Add fast path for minkowski distance with p=1,2

Patch to review:
diff --git a/benchmarks/benchmarks/spatial.py b/benchmarks/benchmarks/spatial.py
index 4fb18ce8a47b..b537a6068b42 100644
--- a/benchmarks/benchmarks/spatial.py
+++ b/benchmarks/benchmarks/spatial.py
@@ -347,24 +347,26 @@ class Xdist(Benchmark):
     def setup(self, num_points, metric):
         np.random.seed(123)
         self.points = np.random.random_sample((num_points, 3))
-        # use an equal weight vector to satisfy those metrics
-        # that require weights
-        if metric == 'wminkowski':
-            self.w = np.ones(3)
+        if metric == 'minkowski':
+            # Default p=2 is just the euclidean metric
+            self.kwargs = {'p': 3.0}
+        elif metric == 'wminkowski':
+            # use an equal weight vector since weights are required
+            self.kwargs = {'w': np.ones(3)}
         else:
-            self.w = None
+            self.kwargs = {}
 
     def time_cdist(self, num_points, metric):
         """Time scipy.spatial.distance.cdist over a range of input data
         sizes and metrics.
         """
-        distance.cdist(self.points, self.points, metric, w=self.w)
+        distance.cdist(self.points, self.points, metric, **self.kwargs)
 
     def time_pdist(self, num_points, metric):
         """Time scipy.spatial.distance.pdist over a range of input data
         sizes and metrics.
         """
-        distance.pdist(self.points, metric, w=self.w)
+        distance.pdist(self.points, metric, **self.kwargs)
 
 
 class ConvexHullBench(Benchmark):
diff --git a/scipy/spatial/src/distance_impl.h b/scipy/spatial/src/distance_impl.h
index 5b2c7154bcd6..314579bbede2 100644
--- a/scipy/spatial/src/distance_impl.h
+++ b/scipy/spatial/src/distance_impl.h
@@ -425,6 +425,78 @@ compute_mean_vector(double *res, const double *X, npy_intp num_rows, const npy_i
 }
 #endif
 
+#define DEFINE_CDIST(name, type) \
+    static int cdist_ ## name ## _ ## type(const type *XA, const type *XB, \
+                                           double *dm,                     \
+                                           const npy_intp num_rowsA,       \
+                                           const npy_intp num_rowsB,       \
+                                           const npy_intp num_cols)        \
+    {                                                                      \
+        Py_ssize_t i, j;                                                   \
+        for (i = 0; i < num_rowsA; ++i) {                                  \
+            const type *u = XA + num_cols * i;                             \
+            for (j = 0; j < num_rowsB; ++j, ++dm) {                        \
+                const type *v = XB + num_cols * j;                         \
+                *dm = name ## _distance_ ## type(u, v, num_cols);          \
+            }                                                              \
+        }                                                                  \
+        return 0;\
+    }
+
+DEFINE_CDIST(bray_curtis, double)
+DEFINE_CDIST(canberra, double)
+DEFINE_CDIST(chebyshev, double)
+DEFINE_CDIST(city_block, double)
+DEFINE_CDIST(euclidean, double)
+DEFINE_CDIST(jaccard, double)
+DEFINE_CDIST(jensenshannon, double)
+DEFINE_CDIST(sqeuclidean, double)
+
+DEFINE_CDIST(dice, char)
+DEFINE_CDIST(jaccard, char)
+DEFINE_CDIST(kulsinski, char)
+DEFINE_CDIST(rogerstanimoto, char)
+DEFINE_CDIST(russellrao, char)
+DEFINE_CDIST(sokalmichener, char)
+DEFINE_CDIST(sokalsneath, char)
+DEFINE_CDIST(yule, char)
+
+
+#define DEFINE_PDIST(name, type) \
+    static int pdist_ ## name ## _ ## type(const type *X, double *dm,       \
+                                           const npy_intp num_rows,         \
+                                           const npy_intp num_cols)         \
+    {                                                                       \
+        Py_ssize_t i, j;                                                    \
+        double *it = dm;                                                    \
+        for (i = 0; i < num_rows; ++i) {                                    \
+            const type *u = X + num_cols * i;                               \
+            for (j = i + 1; j < num_rows; ++j, it++) {                      \
+                const type *v = X + num_cols * j;                           \
+                *it = name ## _distance_ ## type(u, v, num_cols);           \
+            }                                                               \
+        }                                                                   \
+        return 0; \
+    }
+
+DEFINE_PDIST(bray_curtis, double)
+DEFINE_PDIST(canberra, double)
+DEFINE_PDIST(chebyshev, double)
+DEFINE_PDIST(city_block, double)
+DEFINE_PDIST(euclidean, double)
+DEFINE_PDIST(jaccard, double)
+DEFINE_PDIST(jensenshannon, double)
+DEFINE_PDIST(sqeuclidean, double)
+
+DEFINE_PDIST(dice, char)
+DEFINE_PDIST(jaccard, char)
+DEFINE_PDIST(kulsinski, char)
+DEFINE_PDIST(rogerstanimoto, char)
+DEFINE_PDIST(russellrao, char)
+DEFINE_PDIST(sokalmichener, char)
+DEFINE_PDIST(sokalsneath, char)
+DEFINE_PDIST(yule, char)
+
 static NPY_INLINE int
 pdist_mahalanobis(const double *X, double *dm, const npy_intp num_rows,
                   const npy_intp num_cols, const double *covinv)
@@ -499,6 +571,12 @@ pdist_minkowski(const double *X, double *dm, npy_intp num_rows,
                 const npy_intp num_cols, const double p)
 {
     npy_intp i, j;
+    if (p == 1.0) {
+        return pdist_city_block_double(X, dm, num_rows, num_cols);
+    }
+    if (p == 2.0) {
+        return pdist_euclidean_double(X, dm, num_rows, num_cols);
+    }
 
     for (i = 0; i < num_rows; ++i) {
         const double *u = X + (num_cols * i);
@@ -693,6 +771,12 @@ cdist_minkowski(const double *XA, const double *XB, double *dm,
                 const npy_intp num_cols, const double p)
 {
     npy_intp i, j;
+    if (p == 1.0) {
+        return cdist_city_block_double(XA, XB, dm, num_rowsA, num_rowsB, num_cols);
+    }
+    if (p == 2.0) {
+        return cdist_euclidean_double(XA, XB, dm, num_rowsA, num_rowsB, num_cols);
+    }
 
     for (i = 0; i < num_rowsA; ++i) {
         const double *u = XA + (num_cols * i);
@@ -757,75 +841,3 @@ cdist_hamming_char(const char *XA, const char *XB, double *dm,
     }
     return 0;
 }
-
-#define DEFINE_CDIST(name, type) \
-    static int cdist_ ## name ## _ ## type(const type *XA, const type *XB, \
-                                           double *dm,                     \
-                                           const npy_intp num_rowsA,       \
-                                           const npy_intp num_rowsB,       \
-                                           const npy_intp num_cols)        \
-    {                                                                      \
-        Py_ssize_t i, j;                                                   \
-        for (i = 0; i < num_rowsA; ++i) {                                  \
-            const type *u = XA + num_cols * i;                             \
-            for (j = 0; j < num_rowsB; ++j, ++dm) {                        \
-                const type *v = XB + num_cols * j;                         \
-                *dm = name ## _distance_ ## type(u, v, num_cols);          \
-            }                                                              \
-        }                                                                  \
-        return 0;\
-    }
-
-DEFINE_CDIST(bray_curtis, double)
-DEFINE_CDIST(canberra, double)
-DEFINE_CDIST(chebyshev, double)
-DEFINE_CDIST(city_block, double)
-DEFINE_CDIST(euclidean, double)
-DEFINE_CDIST(jaccard, double)
-DEFINE_CDIST(jensenshannon, double)
-DEFINE_CDIST(sqeuclidean, double)
-
-DEFINE_CDIST(dice, char)
-DEFINE_CDIST(jaccard, char)
-DEFINE_CDIST(kulsinski, char)
-DEFINE_CDIST(rogerstanimoto, char)
-DEFINE_CDIST(russellrao, char)
-DEFINE_CDIST(sokalmichener, char)
-DEFINE_CDIST(sokalsneath, char)
-DEFINE_CDIST(yule, char)
-
-
-#define DEFINE_PDIST(name, type) \
-    static int pdist_ ## name ## _ ## type(const type *X, double *dm,       \
-                                           const npy_intp num_rows,         \
-                                           const npy_intp num_cols)         \
-    {                                                                       \
-        Py_ssize_t i, j;                                                    \
-        double *it = dm;                                                    \
-        for (i = 0; i < num_rows; ++i) {                                    \
-            const type *u = X + num_cols * i;                               \
-            for (j = i + 1; j < num_rows; ++j, it++) {                      \
-                const type *v = X + num_cols * j;                           \
-                *it = name ## _distance_ ## type(u, v, num_cols);           \
-            }                                                               \
-        }                                                                   \
-        return 0; \
-    }
-
-DEFINE_PDIST(bray_curtis, double)
-DEFINE_PDIST(canberra, double)
-DEFINE_PDIST(chebyshev, double)
-DEFINE_PDIST(city_block, double)
-DEFINE_PDIST(euclidean, double)
-DEFINE_PDIST(jaccard, double)
-DEFINE_PDIST(jensenshannon, double)
-DEFINE_PDIST(sqeuclidean, double)
-
-DEFINE_PDIST(dice, char)
-DEFINE_PDIST(jaccard, char)
-DEFINE_PDIST(kulsinski, char)
-DEFINE_PDIST(rogerstanimoto, char)
-DEFINE_PDIST(russellrao, char)
-DEFINE_PDIST(sokalmichener, char)
-DEFINE_PDIST(sokalsneath, char)
-DEFINE_PDIST(yule, char)
diff --git a/scipy/spatial/tests/test_distance.py b/scipy/spatial/tests/test_distance.py
index aae361a29705..6cf82f4afba7 100644
--- a/scipy/spatial/tests/test_distance.py
+++ b/scipy/spatial/tests/test_distance.py
@@ -442,28 +442,13 @@ def test_cdist_euclidean_random_unicode(self):
         Y2 = wcdist_no_const(X1, X2, 'test_euclidean')
         _assert_within_tol(Y1, Y2, eps, verbose > 2)
 
-    def test_cdist_minkowski_random_p3d8(self):
+    @pytest.mark.parametrize("p", [1.0, 1.23, 2.0, 3.8, 4.6])
+    def test_cdist_minkowski_random(self, p):
         eps = 1e-07
         X1 = eo['cdist-X1']
         X2 = eo['cdist-X2']
-        Y1 = wcdist_no_const(X1, X2, 'minkowski', p=3.8)
-        Y2 = wcdist_no_const(X1, X2, 'test_minkowski', p=3.8)
-        _assert_within_tol(Y1, Y2, eps, verbose > 2)
-
-    def test_cdist_minkowski_random_p4d6(self):
-        eps = 1e-07
-        X1 = eo['cdist-X1']
-        X2 = eo['cdist-X2']
-        Y1 = wcdist_no_const(X1, X2, 'minkowski', p=4.6)
-        Y2 = wcdist_no_const(X1, X2, 'test_minkowski', p=4.6)
-        _assert_within_tol(Y1, Y2, eps, verbose > 2)
-
-    def test_cdist_minkowski_random_p1d23(self):
-        eps = 1e-07
-        X1 = eo['cdist-X1']
-        X2 = eo['cdist-X2']
-        Y1 = wcdist_no_const(X1, X2, 'minkowski', p=1.23)
-        Y2 = wcdist_no_const(X1, X2, 'test_minkowski', p=1.23)
+        Y1 = wcdist_no_const(X1, X2, 'minkowski', p=p)
+        Y2 = wcdist_no_const(X1, X2, 'test_minkowski', p=p)
         _assert_within_tol(Y1, Y2, eps, verbose > 2)
 
     def test_cdist_cosine_random(self):


Code review comments:
diff_hunk: @@ -347,24 +347,26 @@ class Xdist(Benchmark):
     def setup(self, num_points, metric):
         np.random.seed(123)
         self.points = np.random.random_sample((num_points, 3))
-        # use an equal weight vector to satisfy those metrics
-        # that require weights
-        if metric == 'wminkowski':
-            self.w = np.ones(3)
+        if metric == 'minkowski':
+            # Default p=2 is just the euclidean metric
+            self.kwargs = {'p': 3.0}
path: benchmarks/benchmarks/spatial.py
text: @user1:
Changing this default will make it hard to compare historical performance. Maybe we can add a special case 'minkowski-L3' and have that one add the `p=3` kwarg?

@author:
I thought asv always used the current benchmark source when doing bench-compare?

Even if that's not the case, prior to this change, `p=2` and `p=3` are indistinguishable in performance. So, this shouldn't actually change the result of the benchmark.