Repository: scikit-learn/scikit-learn
Pull request: FIX Raise error when n_neighbours >= n_samples / 2  in manifold.trustworthiness
#### Reference Issues/PRs
Resolves #18832 (stalled)
Fixes #18567

#### What does this implement/fix? Explain your changes.
PR #18832: Added warning to `manifold.trustworthiness` when `n_neighbors > n_features`.
This PR: Improved tests and wording and addressed reviewer comments. 

#### Any other comments?

Problem statement:
bug: sklearn.manifold.turstworthiness outputs values larger than 1
The trustworthiness computation might be erroneous as for some cases I can see that it is providing values that are larger than 1.

Here is a code snippet that produces value larger than 1 in my machine:

```
import numpy as np
from sklearn.manifold import trustworthiness

#import sklearn; sklearn.show_versions()

np.random.seed(5000)

X_train = np.random.rand(7,4)
Y_train = np.random.rand(7,2)

tt = trustworthiness(X_train, Y_train, n_neighbors=5)

print('Computed Trustworthiness:', tt)
```
The output is:
```
Computed Trustworthiness: 1.0857142857142856
```


System:
    python: 3.6.9 (default, Jul 17 2020, 12:50:27)  [GCC 8.4.0]
executable: /usr/bin/python3
   machine: Linux-5.4.0-47-generic-x86_64-with-Ubuntu-18.04-bionic

Python dependencies:
       pip: 20.2.3
setuptools: 41.4.0
   sklearn: 0.22.2.post1
     numpy: 1.17.4
     scipy: 1.4.1
    Cython: None
    pandas: 1.0.5
matplotlib: 3.3.0
    joblib: 0.11

Built with OpenMP: True

Hints from issue discussion:
I get the same result when I run it. I wonder if the implementation of the algorithm is correct or not.

The way it's written in the code doesn't exactly match the description in the paper. Specifically, the code has a `max` function applied where the paper doesn't. But simply removing it didn't fix the problem (the result became negative).

The paper specifies that the inner summation should be over the data vectors that are among the `k` data vectors closest to `x_i` in the original space but not among the `k` data vectors closest to `x_i` after projection. I'm having trouble understanding the middle section of the code, so it's unclear to me whether the code is correctly selecting these data vectors.

Here is a link to the relevant paper: [link](https://www.researchgate.net/publication/2365851_Neighborhood_Preservation_in_Nonlinear_Projection_Methods_An_Experimental_Study).

The implementation of the function is pasted below for reference:

    def trustworthiness(X, X_embedded, *, n_neighbors=5, metric='euclidean'):
        dist_X = pairwise_distances(X, metric=metric)
        if metric == 'precomputed':
            dist_X = dist_X.copy()
        # we set the diagonal to np.inf to exclude the points themselves from
        # their own neighborhood
        np.fill_diagonal(dist_X, np.inf)
        ind_X = np.argsort(dist_X, axis=1)
        # `ind_X[i]` is the index of sorted distances between i and other samples
        ind_X_embedded = NearestNeighbors(n_neighbors=n_neighbors).fit(
                X_embedded).kneighbors(return_distance=False)

        # We build an inverted index of neighbors in the input space: For sample i,
        # we define `inverted_index[i]` as the inverted index of sorted distances:
        # inverted_index[i][ind_X[i]] = np.arange(1, n_sample + 1)
        n_samples = X.shape[0]
        inverted_index = np.zeros((n_samples, n_samples), dtype=int)
        ordered_indices = np.arange(n_samples + 1)
        inverted_index[ordered_indices[:-1, np.newaxis],
                       ind_X] = ordered_indices[1:]
        ranks = inverted_index[ordered_indices[:-1, np.newaxis],
                               ind_X_embedded] - n_neighbors
        t = np.sum(ranks[ranks > 0])
        t = 1.0 - t * (2.0 / (n_samples * n_neighbors *
                              (2.0 * n_samples - 3.0 * n_neighbors - 1.0)))
        return t
I think I found the problem. There is no bug in the code. The bug is in the metric itself.

The metric requires `(2.0 * n_samples - 3.0 * n_neighbors - 1.0) > 0`, which translates to `n_neighbors < (2.0 * n_samples - 1.0) / 3.0`.
For the example above it requires `n_neighbors < (2.0*7-1.0)/3.0 = 4.33`.

In the paper, it is written that "for clarity, we have only included the scaling for neighborhoods of size `k<N/2`". I think we should add this line to the doc.

Or a check in the code to ensure that `k<N/2`. The reference code above can be modified as:

```
def trustworthiness(X, X_embedded, *, n_neighbors=5, metric='euclidean'):
    if n_neighbors >= X.shape[0]/2:
        warning("n_neighbors should be less than' + str(X.shape[0]/2) + 'given '+ str(n_neighbors))
    #other part of the codes
    return t
```

> The way it's written in the code doesn't exactly match the description in the paper. Specifically, the code has a max function applied where the paper doesn't. But simply removing it didn't fix the problem (the result became negative).

The max function is implicit in the paper. It is defined as an intersection of two sets in original and projected data which translates to the max function in practice.
Yes, please.

And perhaps a warning if the algorithm is run on a dataset with such large
n_neighbors!

I will take this. Does just the documentation need changing or does the suggested python change need to also be included?
I think I found the bug and I believe it is unrelated to the value of `k` being greater than `N/2` which shouldn't cause any issue. The problem actually comes from the following line (551):
`inverted_index = np.zeros((n_samples, n_samples), dtype=int)`
Here the array is initialized with `dtype=int`. Further down in the code a sum is performed using `np.sum`. In Numpy's docs it is mentionned that `np.sum` returnes a result with the same `dtype` as the provided input, meaning that in our case the result of the sum is an `int32` with a max value of 2147483648 which is easily reached in some cases. If the result of the sum exceeds this value, the variable `t` becomes negative. As a consequence. the final result becomes >1 which is erroneous.
A simple fix would be to make `inverted_index` an `int64` array which should leave enough room for the sum to never overflow. An additional check could be added to verify that the value stored in `t` is always positive, but that's not very necessary IMO.
I'm willing to submit a PR if my explanation is convincing enough.

Commit message: Added changelog entry

Patch to review:
diff --git a/doc/whats_new/v1.1.rst b/doc/whats_new/v1.1.rst
index 6c03456bcaf0c..59de4efda428e 100644
--- a/doc/whats_new/v1.1.rst
+++ b/doc/whats_new/v1.1.rst
@@ -284,7 +284,7 @@ Changelog
     deprecated.
   - the default value of the `batch_size` parameter of both will change from 3 to 256
     in version 1.3.
-  
+
   :pr:`18975` by :user:`Jérémie du Boisberranger <jeremiedbb>`.
 
 - |Enhancement| :func:`decomposition.dict_learning`, :func:`decomposition.dict_learning_online`
@@ -713,6 +713,12 @@ Changelog
   :pr:`21805` by :user:`Jeanne Dionisi <jeannedionisi>` and :pr:`22685` by
   :user:`Meekail Zain <micky774>`.
 
+- |Enhancement| Add a warning for :func:`manifold.trustworthiness` if
+  `n_neighbours` is less than `n_samples / 2` to ensure a correct
+  support for the function.
+  :pr:`18832` by :user:`Hong Shao Yang <hongshaoyang>` and :pr:`23033` by
+  :user:`Meekail Zain <micky774>`.
+
 :mod:`sklearn.metrics`
 ......................
 
diff --git a/sklearn/manifold/_t_sne.py b/sklearn/manifold/_t_sne.py
index 54ebc1e66d4ae..421788e1185c0 100644
--- a/sklearn/manifold/_t_sne.py
+++ b/sklearn/manifold/_t_sne.py
@@ -474,7 +474,9 @@ def trustworthiness(X, X_embedded, *, n_neighbors=5, metric="euclidean"):
         Embedding of the training data in low-dimensional space.
 
     n_neighbors : int, default=5
-        Number of neighbors k that will be considered.
+        The number of neighbors that will be considered. Should be fewer than
+        `n_samples / 2` to ensure the trustworthiness to lies within [0, 1], as
+        mentioned in [1]_. A warning will be raised otherwise.
 
     metric : str or callable, default='euclidean'
         Which metric to use for computing pairwise distances between samples
@@ -491,7 +493,23 @@ def trustworthiness(X, X_embedded, *, n_neighbors=5, metric="euclidean"):
     -------
     trustworthiness : float
         Trustworthiness of the low-dimensional embedding.
+
+    References
+    ----------
+    .. [1] Jarkko Venna and Samuel Kaski. 2001. Neighborhood
+           Preservation in Nonlinear Projection Methods: An Experimental Study.
+           In Proceedings of the International Conference on Artificial Neural Networks
+           (ICANN '01). Springer-Verlag, Berlin, Heidelberg, 485-491.
+
+    .. [2] Laurens van der Maaten. Learning a Parametric Embedding by Preserving
+           Local Structure. Proceedings of the Twelth International Conference on
+           Artificial Intelligence and Statistics, PMLR 5:384-391, 2009.
     """
+    n_samples = X.shape[0]
+    if n_neighbors >= n_samples / 2:
+        warnings.warn(
+            f"n_neighbors ({n_neighbors}) should be less than {n_samples / 2}"
+        )
     dist_X = pairwise_distances(X, metric=metric)
     if metric == "precomputed":
         dist_X = dist_X.copy()
@@ -509,7 +527,6 @@ def trustworthiness(X, X_embedded, *, n_neighbors=5, metric="euclidean"):
     # We build an inverted index of neighbors in the input space: For sample i,
     # we define `inverted_index[i]` as the inverted index of sorted distances:
     # inverted_index[i][ind_X[i]] = np.arange(1, n_sample + 1)
-    n_samples = X.shape[0]
     inverted_index = np.zeros((n_samples, n_samples), dtype=int)
     ordered_indices = np.arange(n_samples + 1)
     inverted_index[ordered_indices[:-1, np.newaxis], ind_X] = ordered_indices[1:]
diff --git a/sklearn/manifold/tests/test_t_sne.py b/sklearn/manifold/tests/test_t_sne.py
index 0fe70d6dd372d..2f63b241e676d 100644
--- a/sklearn/manifold/tests/test_t_sne.py
+++ b/sklearn/manifold/tests/test_t_sne.py
@@ -1,3 +1,4 @@
+import numbers
 import sys
 from io import StringIO
 import numpy as np
@@ -33,6 +34,7 @@
 from sklearn.metrics.pairwise import pairwise_distances
 from sklearn.metrics.pairwise import manhattan_distances
 from sklearn.metrics.pairwise import cosine_distances
+from sklearn.utils.validation import check_scalar
 
 
 x = np.linspace(0, 1, 10)
@@ -277,6 +279,22 @@ def test_trustworthiness():
     assert_almost_equal(trustworthiness(X, X_embedded, n_neighbors=1), 0.2)
 
 
+def test_trustworthiness_warning():
+    regex = "n_neighbors .+ should be less than .+"
+    X = np.random.rand(7, 4)
+    X_embedded = np.random.rand(7, 2)
+    with pytest.warns(UserWarning, match=regex):
+        trustworthiness(X, X_embedded, n_neighbors=5)
+
+    check_scalar(
+        trustworthiness(X, X_embedded, n_neighbors=3),
+        "trustworthiness",
+        numbers.Real,
+        min_val=0,
+        max_val=1,
+    )
+
+
 @pytest.mark.parametrize("method", ["exact", "barnes_hut"])
 @pytest.mark.parametrize("init", ("random", "pca"))
 def test_preserve_trustworthiness_approximately(method, init):


Code review comments:
diff_hunk: @@ -491,7 +493,23 @@ def trustworthiness(X, X_embedded, *, n_neighbors=5, metric="euclidean"):
     -------
     trustworthiness : float
         Trustworthiness of the low-dimensional embedding.
+
+    References
+    ----------
+    .. [1] Jarkko Venna and Samuel Kaski. 2001. Neighborhood
+           Preservation in Nonlinear Projection Methods: An Experimental Study.
+           In Proceedings of the International Conference on Artificial Neural Networks
+           (ICANN '01). Springer-Verlag, Berlin, Heidelberg, 485-491.
+
+    .. [2] Laurens van der Maaten. Learning a Parametric Embedding by Preserving
+           Local Structure. Proceedings of the Twelth International Conference on
+           Artificial Intelligence and Statistics, PMLR 5:384-391, 2009.
     """
+    n_samples = X.shape[0]
+    if n_neighbors >= n_samples / 2:
+        warnings.warn(
+            f"n_neighbors ({n_neighbors}) should be less than {n_samples / 2}"
path: sklearn/manifold/_t_sne.py
text: @user1:
```suggestion
            f"n_neighbors ({n_neighbors}) should be less than n_samples / 2 ({n_samples / 2})"
```

diff_hunk: @@ -277,6 +279,22 @@ def test_trustworthiness():
     assert_almost_equal(trustworthiness(X, X_embedded, n_neighbors=1), 0.2)
 
 
+def test_trustworthiness_warning():
+    regex = "n_neighbors .+ should be less than .+"
+    X = np.random.rand(7, 4)
+    X_embedded = np.random.rand(7, 2)
+    with pytest.warns(UserWarning, match=regex):
+        trustworthiness(X, X_embedded, n_neighbors=5)
+
+    check_scalar(
+        trustworthiness(X, X_embedded, n_neighbors=3),
path: sklearn/manifold/tests/test_t_sne.py
text: @user1:
I think it's better to assert this directly, rather than using `check_scalar`.