Scipy clustering; use the physics Minkowski metric?

So this morning I learnt that the Minkowski metric does not always mean;

See wolfram for details.

Apparently in scipy it is just a p-norm. Scipy has an option to weight the p-norm, but only with positive weights, so that cannot achieve the relativistic Minkowski metric.

I would like to do hierarchical clustering on points in relativistic 4 dimensional space. For two points;

a = [a_time, a_x, a_y, a_z]

b = [b_time, b_x, b_y, b_z]

The distance between them should be;

invarient_s(a, b) = sqrt(-(a_time-b_time)^2 + (a_x-b_x)^2 + (a_y-b_y)^2 + (a_z-b_z)^2)

I'm working in python and ideally using scipy's fcluster. Before I go and write my own clustering is there anyway to get this metric in fcluster? Can I add to the list of available metrics?

Edit; it appears only fclusterdata supports metrics in the first place.

The bad news is that indeed built-in metrics (and especially the one named Minkowski) don't support negative weights. I suspect the reason for this is that in a proper metric you can only have d(x,y) = 0 if and only if x = y, which is violated by the Minkowski metric. This is probably the reason for the lack of support for negative weights in any of the weighted metrics in scipy, see also remarks in this github thread.

The good news is that the documentation of scipy.cluster.hierarchy.fclusterdata is buggy (now fixed in master), because it claimed

metric: str, optional

    The distance metric for calculating pairwise distances.
    See distance.pdist for descriptions and linkage to verify
    compatibility with the linkage method.

Whereas the actual implementation of fclusterdata simply passes the metric input parameter along to pdist, which allows custom callables to be passed as metric:

metric: str or function, optional

Sure enough, we can define our own Minkowski metric function and pass that on to fclusterdata, but we have to make sure that all the points are spatially separated, otherwise we get complex distances and pdist will loudly fail (complaining about "finite" data, because np.sqrt when given a negative number will return nan, and nan fails the np.isfinite check in linkage). With this reasonable caveat something like the following works:

from scipy.cluster.hierarchy import fclusterdata 
from numpy.random import default_rng  # only for dummy data 
 
# generate random data, use new random machinery for best practices 
N = 10 
rng = default_rng() 
X = rng.random((N, 4)) * [0.01, 1, 1, 1]  # make them all space-like 
 
def physical_minkowski(v1, v2): 
    """Return the proper Minkowski metric for 4-vectors with signature -+++"""
    return np.sqrt(([-1, 1, 1, 1] * v1).dot(v2)) 
 
fclusterdata(X, t=1, metric=physical_minkowski)                               
# returns uninteresting array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=int32)

Since the above function might get called a lot of times, it might make sense to compile it with numba.njit for improved performance. It only needs a small change to make that possible:

import numba

@numba.njit 
def jitted_minkowski(v1, v2): 
    return np.sqrt((np.array([-1, 1, 1, 1]) * v1).dot(v2)) 

I timed both of the above metric functions using IPython's built-in %timeit magic with N = 1000 for a reasonable comparison:

>>> %timeit scipy.spatial.distance.pdist(X, metric=physical_minkowski)
... %timeit scipy.spatial.distance.pdist(X, metric=jitted_minkowski)
2.2 s ± 90.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
385 ms ± 12.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

This means that for larger sets of 4-vectors the JIT-compiled version is 5 times faster, and compilation only has to be done once (you can even cache the compiled function on disk so that you don't have to compile it each time you run your script).