scikit-learn
diff --git a/‎sklearn/metrics/cluster/unsupervised.py
Lines changed: 9 additions & 24 deletions b/‎sklearn/metrics/cluster/unsupervised.py
Lines changed: 9 additions & 24 deletions
@@ -5,19 +5,16 @@
 #          Thierry Guillemot <thierry.guillemot.work@gmail.com>
 # License: BSD 3 clause
 
-from itertools import chain
-from itertools import permutations
-
 import numpy as np
 
-from ...externals.six.moves import xrange
 from ...utils import check_array
 from ...utils import check_consistent_length
 from ...utils import check_random_state
 from ...utils import check_X_y
 from ...utils.fixes import bincount
 from ..pairwise import pairwise_distances
 from ...preprocessing import LabelEncoder
+from .supervised import contingency_matrix
 
 
 def check_number_of_labels(n_labels, n_samples):
@@ -304,27 +301,15 @@ def prediction_strength_score(labels_train, labels_test):
     labels_test = check_array(labels_test, dtype=np.int32, ensure_2d=False,
                               warn_on_dtype=True)
 
-    clusters = set(chain(labels_train, labels_test))
-    n_clusters = len(clusters)
+    n_clusters = max(np.unique(labels_train).shape[0],
+                     np.unique(labels_test).shape[0])
     if n_clusters == 1:
         return 1.0  # by definition
 
-    strength = 1.0
-    for k in clusters:
-        # samples assigned to k-th cluster based on test data
-        samples_test_k = np.flatnonzero(labels_test == k)
-        cluster_test_size = samples_test_k.shape[0]
-
-        if cluster_test_size < 2:
-            continue
-
-        matches = 0
-        for i, j in permutations(xrange(cluster_test_size), 2):
-            ki, kj = samples_test_k[j], samples_test_k[i]
-            if labels_train[ki] == labels_train[kj]:
-                matches += 1
-
-        strength = min(strength, matches / (cluster_test_size *
-                                            (cluster_test_size - 1.)))
+    C = contingency_matrix(labels_train, labels_test)
+    pairs_matching = (C * (C - 1) / 2).sum(axis=0)
+    M = C.sum(axis=0)
+    pairs_total = (M * (M - 1) / 2)
+    nz = pairs_total.nonzero()[0]
 
-    return strength
+    return (pairs_matching[nz] / pairs_total[nz]).min()