PR on scipy (it has been closed)

I should have covered all the problematic pairwise_distances.

Hi,
Thanks for taking a look at this. I read through this PR quickly, and have one major concern: the nested for loops over all the data within the _scipy_wrapper_binary_distances() function. That will be far slower than the current implementation.

SciPy already handles these errors in its C-loop, replacing the values with NaN. It would be much faster and cleaner to do e.g. K[np.isnan(K)] = 0.0 or K[K < 0] = 0, depending on the case.

As far as it's possible, I don't think we should be adding code which duplicates that in scipy.

I agree. But the problem I see with such a post-process check is that we do not have a way to distinguish between constant inputs and genuine bad inputs. For example:

import numpy as np
from sklearn.metrics import pairwise_distances

X1 = np.array([[0,1,0], [0,1,1]])
X2 = np.array([[1,1,1],[0,0,1]])
X3 = X1 + 0.05 # floats, but binary string expected

pairwise_distances(X1,X2, 'yule')
array([[ nan,   2.],
       [ nan,   0.]])
pairwise_distances(X1,X3, 'yule')
array([[ nan,  nan],
       [ nan,  nan]])

Ah, I see. In that case, I think we should wrap the scipy functionality to do both an input validation and output adjustment. Perhaps via function like validate_pairwise_input(X, Y, metric) and filter_pairwise_output(D, metric) or something similar, called at the beginning and end of pairwise_distances.

I read through this PR quickly, and have one major concern: the nested for loops over all the data within the _scipy_wrapper_binary_distances() function. That will be far slower than the current implementation.

@giorgiop can you please run a quick benchmark to check the impact of your current solution on the runtime performance?

Here you are.

import numpy as np
from scipy.spatial.distance import pdist
from sklearn.metrics.pairwise import pairwise_distances

X = np.random.RandomState(0).randint(0, 2, size=(100, 500))

%timeit -n 10 pdist(X, metric='yule')
10 loops, best of 3: 21.2 ms per loop

# Master branch
%timeit -n 10 pairwise_distances(X, metric='yule')
10 loops, best of 3: 21.4 ms per loop

# This branch
%timeit -n 10 pairwise_distances(X, metric='yule')
10 loops, best of 3: 496 ms per loop

I think @jakevdp's solution is the way to go.

Ah, I see. In that case, I think we should wrap the scipy functionality to do both an input validation and output adjustment. Perhaps via function like validate_pairwise_input(X, Y, metric) and filter_pairwise_output(D, metric) or something similar, called at the beginning and end of pairwise_distances.

Unless we call scipy's distance on each pair of vectors (current slow solution), we cannot catch the errors. But we can validate the input, that is, check if the vectors are constant, and then fix the output accordingly.

I think @jakevdp's solution is the way to go.

Indeed.

But we can validate the input, that is, check if the vectors are constant, and then fix the output accordingly.

That's worth trying and benchmarking to see what is the resulting overhead of this strategy.

New benchmark. Slow down only x2

# This branch
%timeit -n 50 pairwise_distances(X, metric='yule')
50 loops, best of 3: 43.4 ms per loop

@jakevdp auxiliary function is avoided, thanks.

# This branch
%timeit -n 50 pairwise_distances(X, metric='yule')
50 loops, best of 3: 42.4 ms per loop

I'm still a bit worried about this conceptually: that is, if the stated definition of, say, yule distance evaluates to 0/0 for a given input, should we really set the NaN result to zero in all cases? I'm not convinced that this choice won't cause problems in other applications. Thoughts, @ogrisel?

For yule, you only get 0/0 if both vectors are constant False, right? [there is not wikipedia entry for this distance?]
Should we warn on constant rows?

For yule, you only get 0/0 if both vectors are constant False, right? [there is not wikipedia entry for this distance?]

0.16 doc (it was wrong in 0.14, still the first result from Google). I think we are dividing by 0 anytime the two vectors are both constant one vector is constant. But let me test (on master):

>>> import numpy as np
>>> from scipy.spatial.distance import cdist
>>> from sklearn.metrics.pairwise import pairwise_distances

>>> X1 = np.zeros((3, 3))
>>> X1[1, :] = 1
>>> X1
array([[ 0.,  0.,  0.],
       [ 1.,  1.,  1.],
       [ 0.,  0.,  0.]])
>>> X2 = np.diag([1, 1, 1])

>>> cdist(X1, X1, 'yule')
array([[ nan,  nan,  nan],
       [ nan,  nan,  nan],
       [ nan,  nan,  nan]])
>>> cdist(X1, X2, 'yule')
array([[ nan,  nan,  nan],
       [ nan,  nan,  nan],
       [ nan,  nan,  nan]])
>>> pairwise_distances(X1, X1, metric='yule')
array([[ 0,  nan,  nan],
       [ nan,  0,  nan],
       [ nan,  nan,  0]])
>>> pairwise_distances(X1, X2, metric='yule')
array([[ nan,  nan,  nan],
       [ nan,  nan,  nan],
       [ nan,  nan,  nan]])

Summary of the changes.

• preprocessig.data._handle_zeros_in_scale now handles also "almost" constant vectors, the ones for which scale is very close to 0 in term of machine precision. This extend the solution for numerical issues while scaling.
• proprocessing.data.normalize can disable this feature by a new param. This is necessary in computing _correlation_distances: we do not want to handle silently those numerical issues all over the code base, but only when required by the application.
• the helper _clip_to_zero set to 0 distances when negative or positive but close to machine precision. If requested, the function also set to 0 the diagonal when X is Y (distance between the vector and itself)
• _correlation_similarity is now a built-in. It's all the same of scipy.spatial but it calls _clip_to_zero internally to cut negative and small positive distances. It still return nan when called on constant vectors (unless they are the same vector), accordingly to scipy semantic. When called with support_constants=True, constant vectors are taken into account. Remark: I did not implement an agnostic filter for nan values: garbage in, garbage out is still in place. We only want to support constant input (for which the distance would be undefined) and numerical issues due to floating point arithmetic.

With regard to binary distances, bad output is filtered only when the function is not well-defined on the input vectors:

• 'yule' => at least one vector is constant
• 'dice' and 'sokalsneath' => both vectors are 0

Could you please have a look and tell me what you think? @jakevdp

The PR should address also #5772 (comment)

The touch to test_gmm apparently fixed a non-deterministic test that was breaking under python=2.6, which I am not able to reproduce locally. It appears to me unrelated to this PR. Honestly I am not sure or can say with certainty that this change is the right one, but I haven't experienced problem with that test anymore.

@jakevdp could you have a look at this when you have time? I haven't touched the PR in a while and it's good to see with a fresh mind whether those changes still make sense.

@giorgiop are you closing this just because it's stagnated? or some other reason?

I have never been fully convinced this was the right approach to solve those numerical issues. I believe this fix was judged as required for the last release but that it was left out for lack of time in reviewing. Interest on this seems to have faded out anyway so I wanted to close. I believe a fresh brand new approach may be better than trying to interpret what I did here.

Let me know if you want to me re-open anyway.