Compute sampled KPCA training scores on every input row

pyod/models/kpca.py

@@ -172,6 +172,8 @@ class KPCA(BaseDetector):
     sampling : bool, optional (default=False)
         If True, sampling subset from the dataset is performed only once,
         in order to reduce time complexity while keeping detection performance.
+        The kernel PCA model is fitted on the subset, but training scores and
+        the outlier threshold are computed from all training samples.
 
     subset_size : float in (0., 1.0) or int (0, n_samples), optional (default=20)
         If sampling is True, the size of subset is specified.
@@ -286,6 +288,7 @@ def fit(self, X, y=None):
         self._set_n_classes(y)
 
         # perform subsampling to reduce time complexity
+        X_fit = X
         if self.sampling is True:
             subset_size = self._check_subset_size(X)
             random_state = check_random_state(self.random_state)
@@ -294,17 +297,17 @@ def fit(self, X, y=None):
                 size=subset_size,
                 replace=False,
             )
-            X = X[random_indices, :]
+            X_fit = X[random_indices, :]
 
         # copy the attributes from the sklearn Kernel PCA object
         if self.n_components is None:
-            n_components = X.shape[0]  # use all dimensions
+            n_components = X_fit.shape[0]  # use all dimensions
         else:
             if self.n_components < 1:
                 raise ValueError(
                     f"`n_components` should be >= 1, got: {self.n_components}"
                 )
-            n_components = min(X.shape[0], self.n_components)
+            n_components = min(X_fit.shape[0], self.n_components)
 
         # validate the number of components to be used for outlier detection
         if self.n_selected_components is None:
@@ -337,18 +340,24 @@ def fit(self, X, y=None):
             n_jobs=self.n_jobs,
             random_state=self.random_state,
         )
-        x_transformed = self.kpca.fit_transform(X)
+        x_transformed = self.kpca.fit_transform(X_fit)
+        if self.sampling is True:
+            x_transformed = self.kpca.transform(X)
         x_transformed = x_transformed[:, : self.n_selected_components_]
 
         centerer = self.kpca.get_centerer
         kernel = self.kpca.get_kernel
+        if self.sampling is True:
+            gram_fit_rows = np.mean(kernel(X, X_fit), axis=1)
+        else:
+            gram_fit_rows = centerer.K_fit_rows_
 
         potential = []
         for i in range(X.shape[0]):
             sample = X[i, :].reshape(1, -1)
             potential.append(kernel(sample))
         potential = np.array(potential).squeeze()
-        potential = potential - 2 * centerer.K_fit_rows_ + centerer.K_fit_all_
+        potential = potential - 2 * gram_fit_rows + centerer.K_fit_all_
 
         # reconstruction error
         self.decision_scores_ = potential - np.sum(np.square(x_transformed),

pyod/test/test_kpca.py

@@ -214,6 +214,56 @@ def tearDown(self):
         pass
 
 
+class TestKPCASamplingScores(unittest.TestCase):
+    def test_scores_and_threshold_use_all_training_rows(self):
+        X = np.random.RandomState(7).normal(size=(80, 4))
+        X[-8:] += 3
+        for kernel in ("linear", "rbf"):
+            for subset_size in (12, 0.25, len(X)):
+                with self.subTest(kernel=kernel, subset_size=subset_size):
+                    size = (int(len(X) * subset_size)
+                            if isinstance(subset_size, float)
+                            else subset_size)
+                    indices = np.random.RandomState(42).choice(
+                        len(X), size=size, replace=False)
+                    basis = X[indices]
+                    if kernel == "linear":
+                        K = basis @ basis.T
+                        cross = X @ basis.T
+                        diagonal = np.sum(X ** 2, axis=1)
+                    else:
+                        K = np.exp(-0.25 * np.sum(
+                            (basis[:, None] - basis[None, :]) ** 2, axis=2))
+                        cross = np.exp(-0.25 * np.sum(
+                            (X[:, None] - basis[None, :]) ** 2, axis=2))
+                        diagonal = np.ones(len(X))
+
+                    # Independent centered-kernel eigendecomposition.
+                    row_mean = K.mean(axis=0)
+                    K_centered = (K - row_mean[None, :]
+                                  - row_mean[:, None] + K.mean())
+                    eigenvalues, eigenvectors = np.linalg.eigh(K_centered)
+                    projection = ((cross - cross.mean(axis=1, keepdims=True)
+                                   - row_mean + K.mean())
+                                  @ eigenvectors[:, -2:]
+                                  / np.sqrt(eigenvalues[-2:]))
+                    expected = (diagonal - 2 * cross.mean(axis=1) + K.mean()
+                                - np.sum(projection ** 2, axis=1))
+
+                    clf = KPCA(sampling=True, subset_size=subset_size,
+                               n_components=3, n_selected_components=2,
+                               kernel=kernel, gamma=0.25, eigen_solver="dense",
+                               contamination=0.1, random_state=42).fit(X)
+                    assert_equal(clf.kpca.X_fit_, basis)
+                    assert_allclose(clf.decision_scores_, expected, atol=1e-12)
+                    assert_allclose(clf.decision_scores_,
+                                    clf.decision_function(X), atol=1e-12)
+                    assert_allclose(clf.threshold_,
+                                    np.percentile(expected, 90), atol=1e-12)
+                    assert_equal(clf.labels_, expected > clf.threshold_)
+                    assert_equal(clf.labels_, clf.predict(X))
+
+
 class TestKPCAComponentsBound(unittest.TestCase):
     def setUp(self):
         self.n_train = 200