Repository: scikit-learn/scikit-learn
Pull request: [MRG] ENH Add parameter to sklearn.manifold.TSNE to transition to squaring non-Euclidean distance metrics by default (#12401)
#### Reference Issues/PRs

Fixes #12401
Supersedes and closes #12419

The plan for deprecating legacy squaring behaviour is summarized in [this comment.](https://github.com/scikit-learn/scikit-learn/issues/12401#issuecomment-652893433)

#### What does this implement/fix? Explain your changes.

 * Adds a `square_distances` parameter to TSNE. (`default='legacy'`)
     * `True` enables squaring for all distances metrics
     * `'legacy'` preserves existing behaviour, where distances will only be squared when `metric='euclidean'`. Raises a `FutureWarning` when a non-Euclidean metric is provided to alert the deprecation of legacy behavior. Setting `True` suppresses this warning.
 * Adds a parameter description to the docstring of TSNE
 * Adds tests to check `square_distances` and check that `FutureWarning` is raised properly.
 * Modifies existing tests to include `square_distances=True` to suppress FutureWarnings.

#### Any other comments?

* ~~I was uncertain about the usage of `pytest.mark.parametrize` in my tests. Some existing tests use it, while others iterate through lists of parameters. Style suggestions are appreciated.~~ Answered in feedback.
* ~~I was uncertain about the descriptions for my tests. Some use docstrings while others use a simple comment. Style suggestions are appreciated.~~ Answered in feedback.
* I was uncertain about proposing a timeline for changing the default value and/or deprecating `'legacy'` behaviour. Discussion in the issue (#12401) has become outdated, so I chose version numbers based on the deprecation section of the Contributing docs.

Problem statement:
tSNE with non-Euclidean metric should be able to square the distances
Our tSNE implementation uses squared Euclidean distances by default, but does not square the distances when other metrics, or precomputed data, are provided. We had no certainty about whether the theory underlying tSNE was even valid for non-Euclidean metrics (which was raised in the partial rewrite before version 0.19) but if it is valid I would think using a matrix of squared pairwise distances would be more appropriate than raw distance matrices. I think we should at least facilitate squaring the distance matrix, if not deprecate the default behaviour.

I propose:
* Add a parameter square_distance, setting is default value to 'warn' (or 'legacy')
* Current (legacy) behaviour would square only if metric='euclidean'
* When True, all distance matrices will be squared explicitly, including precomputed
* When False, no distance matrices will be squared
* When 'warn', a FutureWarning will be emitted if metric is not Euclidean: from version 0.23 square_distance will default to True

Hints from issue discussion:
I'd like to work on this.
I'm not sure that we have consensus from the core devs that this is the
right solution, but I'd be happy for you to work on it if you're confident
on what to do and how to test it.

I am not sure if tSNE should be used for non-Euclidean distances. Could not find anything
 in the literature. Nevertheless, if there is agreement to the steps you suggested, then you can take a look at my PR.
I have marked it WIP until a common consensus is reached. Will add pytests and more documentation if the approach seems okay.

Thanks!
> I think we should at least facilitate squaring the distance matrix.

I agree.
Hi Joel,

I have made the respective changes as per your recent comments, but I am struggling with unit tests as there are a few unit tests failing in K-means which are completely unrelated to my changes. I have been trying to solve them for last week but couldn't get any luck. 

Could you suggest me what is the required course of action here? Should I just upload my changes and wait for the other unit tests to be fixed?

Thanks,
Anurag
I have no idea why unit tests might be failing in k-means. There are no such failures in master on our continuous integration servers. I do think you should upload your changes so that we can see what fails on Travis and can help you resolve them. (Please discuss further pull request related discussion at #12524, not here.)
Seeing that #12419 is stalled, I'm interested in picking this up.

I see that the `help wanted` tag was removed by @adrinjalali, though. Is this enhancement still wanted? Thanks! :)
Feel free to take it @joshuacwnewton . The label was removed in 2018, and this issue was just never triaged after that :D
Take
I have a naive question here but why would squaring the distance matrix be useful ?
For euclidean it makes sense because we first compute the squared distance matrix and then take the square root.

Commit message: Fix line length linting issues

`square_distance` -> `square_distances` pushed some lines > 79.

Patch to review:
diff --git a/doc/whats_new/v0.24.rst b/doc/whats_new/v0.24.rst
index 895a9863cada8..36761ea8ef212 100644
--- a/doc/whats_new/v0.24.rst
+++ b/doc/whats_new/v0.24.rst
@@ -161,6 +161,13 @@ Changelog
 - |Enhancement| :class:`isotonic.IsotonicRegression` now accepts 2darray with 1 feature as
   input array. :pr:`17379` by :user:`Jiaxiang <fujiaxiang>`.
 
+:mod:`sklearn.manifold`
+.......................
+
+- |Enhancement| Add `square_distance` parameter to :class:`manifold.TSNE`,
+  which provides backward compatibility during deprecation of legacy squaring
+  behavior. :pr:`17662` by :user:`Joshua Newton <joshuacwnewton>`.
+
 :mod:`sklearn.metrics`
 ......................
 
diff --git a/sklearn/manifold/_t_sne.py b/sklearn/manifold/_t_sne.py
index 2e6c64362775b..0557ecebd8f0b 100644
--- a/sklearn/manifold/_t_sne.py
+++ b/sklearn/manifold/_t_sne.py
@@ -8,6 +8,7 @@
 # * Fast Optimization for t-SNE:
 #   https://cseweb.ucsd.edu/~lvdmaaten/workshops/nips2010/papers/vandermaaten.pdf
 
+import warnings
 from time import time
 import numpy as np
 from scipy import linalg
@@ -597,6 +598,19 @@ class TSNE(BaseEstimator):
 
         .. versionadded:: 0.22
 
+    square_distances : {True, 'legacy'}, default='legacy'
+        Whether TSNE should square the distance values. ``'legacy'`` means
+        that distance values are squared only when ``metric="euclidean"``.
+        ``True`` means that distance values are squared for all metrics.
+
+        .. versionadded:: 0.24
+           Added to provide backward compatibility during deprecation of
+           legacy squaring behavior.
+        .. deprecated:: 0.24
+           Legacy squaring behavior was deprecated in 0.24. The ``'legacy'``
+           value will be removed in 0.26, at which point the default value will
+           change to ``True``.
+
     Attributes
     ----------
     embedding_ : array-like, shape (n_samples, n_components)
@@ -643,7 +657,7 @@ def __init__(self, n_components=2, *, perplexity=30.0,
                  n_iter_without_progress=300, min_grad_norm=1e-7,
                  metric="euclidean", init="random", verbose=0,
                  random_state=None, method='barnes_hut', angle=0.5,
-                 n_jobs=None):
+                 n_jobs=None, square_distances='legacy'):
         self.n_components = n_components
         self.perplexity = perplexity
         self.early_exaggeration = early_exaggeration
@@ -658,6 +672,8 @@ def __init__(self, n_components=2, *, perplexity=30.0,
         self.method = method
         self.angle = angle
         self.n_jobs = n_jobs
+        # TODO Revisit deprecation of square_distances for 0.26-0.28 (#12401)
+        self.square_distances = square_distances
 
     def _fit(self, X, skip_num_points=0):
         """Private function to fit the model using X as training data."""
@@ -666,6 +682,15 @@ def _fit(self, X, skip_num_points=0):
             raise ValueError("'method' must be 'barnes_hut' or 'exact'")
         if self.angle < 0.0 or self.angle > 1.0:
             raise ValueError("'angle' must be between 0.0 - 1.0")
+        if self.square_distances not in [True, 'legacy']:
+            raise ValueError("'square_distances' must be True or 'legacy'.")
+        if self.metric != "euclidean" and self.square_distances is not True:
+            warnings.warn(("'square_distances' has been introduced in 0.24. "
+                           "'legacy' is provided as a setting for backward "
+                           "compatibility purposes. However, 'legacy' will "
+                           "be removed in 0.26, at which point all distance "
+                           "metrics will be squared by default."),
+                          FutureWarning)
         if self.method == 'barnes_hut':
             X = self._validate_data(X, accept_sparse=['csr'],
                                     ensure_min_samples=2,
@@ -715,15 +740,23 @@ def _fit(self, X, skip_num_points=0):
                     print("[t-SNE] Computing pairwise distances...")
 
                 if self.metric == "euclidean":
+                    # Euclidean is squared here, rather than using **= 2,
+                    # because euclidean_distances already calculates
+                    # squared distances, and returns np.sqrt(dist) for
+                    # squared=False.
+                    # Also, Euclidean is slower for n_jobs>1, so don't set here
                     distances = pairwise_distances(X, metric=self.metric,
                                                    squared=True)
                 else:
                     distances = pairwise_distances(X, metric=self.metric,
                                                    n_jobs=self.n_jobs)
 
-                if np.any(distances < 0):
-                    raise ValueError("All distances should be positive, the "
-                                     "metric given is not correct")
+            if np.any(distances < 0):
+                raise ValueError("All distances should be positive, the "
+                                 "metric given is not correct")
+
+            if self.metric != "euclidean" and self.square_distances is True:
+                distances **= 2
 
             # compute the joint probability distribution for the input space
             P = _joint_probabilities(distances, self.perplexity, self.verbose)
@@ -765,7 +798,7 @@ def _fit(self, X, skip_num_points=0):
             # Free the memory used by the ball_tree
             del knn
 
-            if self.metric == "euclidean":
+            if self.square_distances is True or self.metric == "euclidean":
                 # knn return the euclidean distance but we need it squared
                 # to be consistent with the 'exact' method. Note that the
                 # the method was derived using the euclidean method as in the
diff --git a/sklearn/manifold/tests/test_t_sne.py b/sklearn/manifold/tests/test_t_sne.py
index 5e38e1afff592..5ee2b10fd05ca 100644
--- a/sklearn/manifold/tests/test_t_sne.py
+++ b/sklearn/manifold/tests/test_t_sne.py
@@ -285,7 +285,8 @@ def test_preserve_trustworthiness_approximately_with_precomputed_distances():
         D = squareform(pdist(X), "sqeuclidean")
         tsne = TSNE(n_components=2, perplexity=2, learning_rate=100.0,
                     early_exaggeration=2.0, metric="precomputed",
-                    random_state=i, verbose=0, n_iter=500)
+                    random_state=i, verbose=0, n_iter=500,
+                    square_distances=True)
         X_embedded = tsne.fit_transform(D)
         t = trustworthiness(D, X_embedded, n_neighbors=1, metric="precomputed")
         assert t > .95
@@ -325,13 +326,14 @@ def test_too_few_iterations():
     ([[0., -1.], [1., 0.]], ".* positive.*"),
 ])
 def test_bad_precomputed_distances(method, D, retype, message_regex):
-    tsne = TSNE(metric="precomputed", method=method)
+    tsne = TSNE(metric="precomputed", method=method,
+                square_distances=True)
     with pytest.raises(ValueError, match=message_regex):
         tsne.fit_transform(retype(D))
 
 
 def test_exact_no_precomputed_sparse():
-    tsne = TSNE(metric='precomputed', method='exact')
+    tsne = TSNE(metric='precomputed', method='exact', square_distances=True)
     with pytest.raises(TypeError, match='sparse'):
         tsne.fit_transform(sp.csr_matrix([[0, 5], [5, 0]]))
 
@@ -340,7 +342,7 @@ def test_high_perplexity_precomputed_sparse_distances():
     # Perplexity should be less than 50
     dist = np.array([[1., 0., 0.], [0., 1., 0.], [1., 0., 0.]])
     bad_dist = sp.csr_matrix(dist)
-    tsne = TSNE(metric="precomputed")
+    tsne = TSNE(metric="precomputed", square_distances=True)
     msg = "3 neighbors per samples are required, but some samples have only 1"
     with pytest.raises(ValueError, match=msg):
         tsne.fit_transform(bad_dist)
@@ -358,7 +360,7 @@ def test_sparse_precomputed_distance():
     assert sp.issparse(D_sparse)
     assert_almost_equal(D_sparse.A, D)
 
-    tsne = TSNE(metric="precomputed", random_state=0)
+    tsne = TSNE(metric="precomputed", random_state=0, square_distances=True)
     Xt_dense = tsne.fit_transform(D)
 
     for fmt in ['csr', 'lil']:
@@ -371,7 +373,8 @@ def test_non_positive_computed_distances():
     def metric(x, y):
         return -1
 
-    tsne = TSNE(metric=metric, method='exact')
+    # Negative computed distances should be caught even if result is squared
+    tsne = TSNE(metric=metric, method='exact', square_distances=True)
     X = np.array([[0.0, 0.0], [1.0, 1.0]])
     with pytest.raises(ValueError, match="All distances .*metric given.*"):
         tsne.fit_transform(X)
@@ -395,17 +398,19 @@ def test_init_ndarray():
 def test_init_ndarray_precomputed():
     # Initialize TSNE with ndarray and metric 'precomputed'
     # Make sure no FutureWarning is thrown from _fit
-    tsne = TSNE(init=np.zeros((100, 2)), metric="precomputed")
+    tsne = TSNE(init=np.zeros((100, 2)), metric="precomputed",
+                square_distances=True)
     tsne.fit(np.zeros((100, 100)))
 
 
 def test_distance_not_available():
     # 'metric' must be valid.
-    tsne = TSNE(metric="not available", method='exact')
+    tsne = TSNE(metric="not available", method='exact', square_distances=True)
     with pytest.raises(ValueError, match="Unknown metric not available.*"):
         tsne.fit_transform(np.array([[0.0], [1.0]]))
 
-    tsne = TSNE(metric="not available", method='barnes_hut')
+    tsne = TSNE(metric="not available", method='barnes_hut',
+                square_distances=True)
     with pytest.raises(ValueError, match="Metric 'not available' not valid.*"):
         tsne.fit_transform(np.array([[0.0], [1.0]]))
 
@@ -417,6 +422,13 @@ def test_method_not_available():
         tsne.fit_transform(np.array([[0.0], [1.0]]))
 
 
+def test_square_distances_not_available():
+    # square_distances must be True or 'legacy'.
+    tsne = TSNE(square_distances="not_available")
+    with pytest.raises(ValueError, match="'square_distances' must be True or"):
+        tsne.fit_transform(np.array([[0.0], [1.0]]))
+
+
 def test_angle_out_of_range_checks():
     # check the angle parameter range
     for angle in [-1, -1e-6, 1 + 1e-6, 2]:
@@ -428,7 +440,7 @@ def test_angle_out_of_range_checks():
 
 def test_pca_initialization_not_compatible_with_precomputed_kernel():
     # Precomputed distance matrices must be square matrices.
-    tsne = TSNE(metric="precomputed", init="pca")
+    tsne = TSNE(metric="precomputed", init="pca", square_distances=True)
     with pytest.raises(ValueError, match="The parameter init=\"pca\" cannot"
                                          " be used with"
                                          " metric=\"precomputed\"."):
@@ -586,7 +598,7 @@ def test_verbose():
 def test_chebyshev_metric():
     # t-SNE should allow metrics that cannot be squared (issue #3526).
     random_state = check_random_state(0)
-    tsne = TSNE(metric="chebyshev")
+    tsne = TSNE(metric="chebyshev", square_distances=True)
     X = random_state.randn(5, 2)
     tsne.fit_transform(X)
 
@@ -872,13 +884,64 @@ def test_tsne_with_different_distance_metrics():
     for metric, dist_func in zip(metrics, dist_funcs):
         X_transformed_tsne = TSNE(
             metric=metric, n_components=n_components_embedding,
-            random_state=0, n_iter=300).fit_transform(X)
+            random_state=0, n_iter=300, square_distances=True).fit_transform(X)
         X_transformed_tsne_precomputed = TSNE(
             metric='precomputed', n_components=n_components_embedding,
-            random_state=0, n_iter=300).fit_transform(dist_func(X))
+            random_state=0, n_iter=300,
+            square_distances=True).fit_transform(dist_func(X))
         assert_array_equal(X_transformed_tsne, X_transformed_tsne_precomputed)
 
 
+@pytest.mark.parametrize('method', ['exact', 'barnes_hut'])
+@pytest.mark.parametrize('metric', ['euclidean', 'manhattan'])
+@pytest.mark.parametrize('square_distances', [True, 'legacy'])
+@ignore_warnings(category=FutureWarning)
+def test_tsne_different_square_distances(method, metric, square_distances):
+    # Make sure that TSNE works for different square_distances settings
+    # FIXME remove test when square_distances=True becomes the default in 0.26
+    random_state = check_random_state(0)
+    n_components_original = 3
+    n_components_embedding = 2
+
+    # Used to create data with structure; this avoids unstable behavior in TSNE
+    X, _ = make_blobs(n_features=n_components_original,
+                      random_state=random_state)
+    X_precomputed = pairwise_distances(X, metric=metric)
+
+    if metric == 'euclidean' and square_distances == 'legacy':
+        X_precomputed **= 2
+
+    X_transformed_tsne = TSNE(
+        metric=metric, n_components=n_components_embedding,
+        square_distances=square_distances, method=method,
+        random_state=0).fit_transform(X)
+    X_transformed_tsne_precomputed = TSNE(
+        metric='precomputed', n_components=n_components_embedding,
+        square_distances=square_distances, method=method,
+        random_state=0).fit_transform(X_precomputed)
+
+    assert_allclose(X_transformed_tsne, X_transformed_tsne_precomputed)
+
+
+@pytest.mark.parametrize('metric', ['euclidean', 'manhattan'])
+@pytest.mark.parametrize('square_distances', [True, 'legacy'])
+def test_tsne_square_distances_futurewarning(metric, square_distances):
+    # Make sure that a FutureWarning is only raised when a non-Euclidean
+    # metric is specified and square_metric is not set to True.
+    random_state = check_random_state(0)
+
+    X = random_state.randn(5, 2)
+    tsne = TSNE(metric=metric, square_distances=square_distances)
+
+    if metric != 'euclidean' and square_distances is not True:
+        with pytest.warns(FutureWarning, match="'square_distances'.*"):
+            tsne.fit_transform(X)
+    else:
+        with pytest.warns(None) as record:
+            tsne.fit_transform(X)
+        assert not record
+
+
 @pytest.mark.parametrize('method', ['exact', 'barnes_hut'])
 def test_tsne_n_jobs(method):
     """Make sure that the n_jobs parameter doesn't impact the output"""
diff --git a/sklearn/neighbors/tests/test_neighbors_pipeline.py b/sklearn/neighbors/tests/test_neighbors_pipeline.py
index 1f41298aa6256..f8f9472bdac48 100644
--- a/sklearn/neighbors/tests/test_neighbors_pipeline.py
+++ b/sklearn/neighbors/tests/test_neighbors_pipeline.py
@@ -127,9 +127,11 @@ def test_tsne():
             KNeighborsTransformer(n_neighbors=n_neighbors, mode='distance',
                                   metric=metric),
             TSNE(metric='precomputed', perplexity=perplexity,
-                 method="barnes_hut", random_state=42, n_iter=n_iter))
+                 method="barnes_hut", random_state=42, n_iter=n_iter,
+                 square_distances=True))
         est_compact = TSNE(metric=metric, perplexity=perplexity, n_iter=n_iter,
-                           method="barnes_hut", random_state=42)
+                           method="barnes_hut", random_state=42,
+                           square_distances=True)
 
         Xt_chain = est_chain.fit_transform(X)
         Xt_compact = est_compact.fit_transform(X)


Code review comments:
diff_hunk: @@ -872,13 +884,64 @@ def test_tsne_with_different_distance_metrics():
     for metric, dist_func in zip(metrics, dist_funcs):
         X_transformed_tsne = TSNE(
             metric=metric, n_components=n_components_embedding,
-            random_state=0, n_iter=300).fit_transform(X)
+            random_state=0, n_iter=300, square_distances=True).fit_transform(X)
         X_transformed_tsne_precomputed = TSNE(
             metric='precomputed', n_components=n_components_embedding,
-            random_state=0, n_iter=300).fit_transform(dist_func(X))
+            random_state=0, n_iter=300,
+            square_distances=True).fit_transform(dist_func(X))
         assert_array_equal(X_transformed_tsne, X_transformed_tsne_precomputed)
 
 
+@pytest.mark.parametrize('method', ['exact', 'barnes_hut'])
+@pytest.mark.parametrize('metric', ['euclidean', 'manhattan'])
+@pytest.mark.parametrize('square_distances', [True, 'legacy'])
+@ignore_warnings(category=FutureWarning)
+def test_tsne_different_square_distances(method, metric, square_distances):
+    # Make sure that TSNE works for different square_distances settings
+    # FIXME remove test when square_distances=True becomes the default in 0.26
+    random_state = check_random_state(0)
+    n_components_original = 3
+    n_components_embedding = 2
+
+    # Used to create data with structure; this avoids unstable behavior in TSNE
+    X, _ = make_blobs(n_features=n_components_original,
+                      random_state=random_state)
+    X_precomputed = pairwise_distances(X, metric=metric)
+
+    if metric == 'euclidean' and square_distances == 'legacy':
+        X_precomputed **= 2
+
+    X_transformed_tsne = TSNE(
+        metric=metric, n_components=n_components_embedding,
+        square_distances=square_distances, method=method,
+        random_state=0).fit_transform(X)
+    X_transformed_tsne_precomputed = TSNE(
+        metric='precomputed', n_components=n_components_embedding,
+        square_distances=square_distances, method=method,
+        random_state=0).fit_transform(X_precomputed)
+
+    assert_allclose(X_transformed_tsne, X_transformed_tsne_precomputed)
+
+
+@pytest.mark.parametrize('metric', ['euclidean', 'manhattan'])
+@pytest.mark.parametrize('square_distances', [True, 'legacy'])
+def test_tsne_square_distances_futurewarning(metric, square_distances):
+    # Make sure that a FutureWarning is only raised when a non-Euclidean
+    # metric is specified and square_metric is not set to True.
path: sklearn/manifold/tests/test_t_sne.py
text: @user1:
thank you for writing helpful test comments! these are important

@author:
Ah, I've just noticed there's a typo in the test comment you quoted (`square_metric`)! I'll fix that now. :)

diff_hunk: @@ -666,6 +682,15 @@ def _fit(self, X, skip_num_points=0):
             raise ValueError("'method' must be 'barnes_hut' or 'exact'")
         if self.angle < 0.0 or self.angle > 1.0:
             raise ValueError("'angle' must be between 0.0 - 1.0")
+        if self.square_distances not in [True, 'legacy']:
+            raise ValueError("'square_distances' must be True or 'legacy'.")
+        if self.metric != "euclidean" and self.square_distances is not True:
+            warnings.warn(("'square_distances' has been introduced in 0.24. "
+                           "'legacy' is provided as a setting for backward "
+                           "compatibility purposes. However, 'legacy' will "
+                           "be removed in 0.26, at which point all distance "
+                           "metrics will be squared by default."),
path: sklearn/manifold/_t_sne.py
text: @user1:
let's add

> Set square_distances=True to avoid this warning. Note that this parameter will be removed in Y.

@author:
Hmm. This is a small thing and probably not worth worrying too much about, but, do you have any thoughts on https://github.com/scikit-learn/scikit-learn/pull/17662#discussion_r450137420?

@user1:
I'd still add it, because otherwise users have no clue what to do unless they go check out the new docs.

@author:
Thanks! I've re-written the message in a way that I hope addresses both of your feedback. :)

diff_hunk: @@ -161,6 +161,13 @@ Changelog
 - |Enhancement| :class:`isotonic.IsotonicRegression` now accepts 2darray with 1 feature as
   input array. :pr:`17379` by :user:`Jiaxiang <fujiaxiang>`.
 
+:mod:`sklearn.manifold`
+.......................
+
+- |Enhancement| Add `square_distance` parameter to :class:`manifold.TSNE`,
path: doc/whats_new/v0.24.rst
text: @user1:
```suggestion
- |Enhancement| Add `square_distances` parameter to :class:`manifold.TSNE`,
```

@author:
Good catch! :)