scikit-learn · dkobak · May 6, 2025 · May 6, 2025 · May 6, 2025 · May 6, 2025
diff --git a/doc/api_reference.py b/doc/api_reference.py
@@ -694,6 +694,7 @@ def _get_submodule(module_name, submodule_name):
                     "Isomap",
                     "LocallyLinearEmbedding",
                     "MDS",
+                    "ClassicalMDS",
                     "SpectralEmbedding",
-                    "Isomap",
-                    "LocallyLinearEmbedding",
-                    "MDS",
-                    "ClassicalMDS",
-                    "SpectralEmbedding",
+                    "ClassicalMDS",
+                    "Isomap",
+                    "LocallyLinearEmbedding",
+                    "MDS",
+                    "SpectralEmbedding",
-                    "Isomap",
-                    "LocallyLinearEmbedding",
-                    "MDS",
-                    "ClassicalMDS",
-                    "SpectralEmbedding",
+                    "ClassicalMDS",
+                    "Isomap",
+                    "LocallyLinearEmbedding",
+                    "MDS",
+                    "SpectralEmbedding",
                     "TSNE",
                     "locally_linear_embedding",

diff --git a/doc/modules/manifold.rst b/doc/modules/manifold.rst
@@ -489,6 +489,15 @@ coordinates :math:`Z` of the embedded points.
     :align: center
     :scale: 60
 
+Apart from that, there is a version called *classical MDS*, also known as
+*principal coordinates analysis (PCoA)* or *Torgerson's scaling*, and implemented
+in the separate :class:`ClassicalMDS` class. Classical  MDS replaces the stress
+loss function with a different loss function called *strain*, which allows
+exact solution in terms of eigendecomposition of the double-centered dissimilarity
+matrix. If the dissimilarity matrix consists of the pairwise Euclidean distances
+between some vectors, then classical MDS is equivalent to PCA applied to this
+set of vectors.
+
 .. rubric:: References
 
 * `"More on Multidimensional Scaling and Unfolding in R: smacof Version 2"

diff --git a/doc/whats_new/upcoming_changes/sklearn.manifold/31322.feature.rst b/doc/whats_new/upcoming_changes/sklearn.manifold/31322.feature.rst
@@ -0,0 +1,3 @@
+- :class:`manifold.ClassicalMDS` was implemented to perform classical MDS
+  (eigendecomposition of the double-centered distance matrix).
+  By :user:`Dmitry Kobak <dkobak>` and :user:`Meekail Zain <Micky774>`
diff --git a/sklearn/manifold/__init__.py b/sklearn/manifold/__init__.py
@@ -3,6 +3,7 @@
 # Authors: The scikit-learn developers
 # SPDX-License-Identifier: BSD-3-Clause
 
+from ._classical_mds import ClassicalMDS
 from ._isomap import Isomap
 from ._locally_linear import LocallyLinearEmbedding, locally_linear_embedding
 from ._mds import MDS, smacof
@@ -12,6 +13,7 @@
 __all__ = [
     "MDS",
     "TSNE",
+    "ClassicalMDS",
     "Isomap",
     "LocallyLinearEmbedding",
     "SpectralEmbedding",

diff --git a/sklearn/manifold/_classical_mds.py b/sklearn/manifold/_classical_mds.py
@@ -0,0 +1,192 @@
+"""
+Classical multi-dimensional scaling (classical MDS).
+"""
+
+# Authors: The scikit-learn developers
+# SPDX-License-Identifier: BSD-3-Clause
+
+from numbers import Integral
+
+import numpy as np
+from scipy import linalg
+
+from ..base import BaseEstimator, _fit_context
+from ..metrics import pairwise_distances
+from ..utils import check_symmetric
+from ..utils._param_validation import Interval
+from ..utils.extmath import svd_flip
+from ..utils.validation import validate_data
+
+
+class ClassicalMDS(BaseEstimator):
+    """Classical multidimensional scaling.
+
+    Read more in the :ref:`User Guide <multidimensional_scaling>`.
+
+    Parameters
+    ----------
+    n_components : int, default=2
+        Number of embedding dimensions.
+
+    metric : str or callable, default='euclidean'
+        Metric to use for dissimilarity computation. Default is "euclidean".
+
+        If metric is a string, it must be one of the options allowed by
+        `scipy.spatial.distance.pdist` for its metric parameter, or a metric
+        listed in :func:`sklearn.metrics.pairwise.distance_metrics`
+
+        If metric is "precomputed", X is assumed to be a distance matrix and
+        must be square during fit.
+
+        If metric is a callable function, it takes two arrays representing 1D
+        vectors as inputs and must return one value indicating the distance
+        between those vectors. This works for Scipy's metrics, but is less
+        efficient than passing the metric name as a string.
+
+    metric_params : dict, default=None
+        Additional keyword arguments for the dissimilarity computation.
+
+    Attributes
+    ----------
+    embedding_ : ndarray of shape (n_samples, n_components)
+        Stores the position of the dataset in the embedding space.
+
+    dissimilarity_matrix_ : ndarray of shape (n_samples, n_samples)
+        Pairwise dissimilarities between the points.
+
+    eigenvalues_ : ndarray of shape (n_components,)
+        Eigenvalues of the double-centered dissimilarity matrix, corresponding
+        to each of the selected components. They are equal to the squared 2-norms
+        of the `n_components` variables in the embedding space.
+
+    n_features_in_ : int
+        Number of features seen during :term:`fit`.
+
+    feature_names_in_ : ndarray of shape (`n_features_in_`,)
+        Names of features seen during :term:`fit`. Defined only when `X`
+        has feature names that are all strings.
+
+    See Also
+    --------
+    sklearn.decomposition.PCA : Principal component analysis.
+    MDS : Metric and non-metric MDS.
+
+    References
+    ----------
+    .. [1] "Modern Multidimensional Scaling - Theory and Applications" Borg, I.;
+       Groenen P. Springer Series in Statistics (1997)
+
+    Examples
+    --------
+    >>> from sklearn.datasets import load_digits
+    >>> from sklearn.manifold import ClassicalMDS
+    >>> X, _ = load_digits(return_X_y=True)
+    >>> X.shape
+    (1797, 64)
+    >>> cmds = ClassicalMDS(n_components=2)
+    >>> X_emb = cmds.fit_transform(X[:100])
+    >>> X_emb.shape
+    (100, 2)
+    """
+
+    _parameter_constraints: dict = {
+        "n_components": [Interval(Integral, 1, None, closed="left")],
+        "metric": [str, callable],
+        "metric_params": [dict, None],
+    }
+
+    def __init__(
+        self,
+        n_components=2,
+        *,
+        metric="euclidean",
+        metric_params=None,
+    ):
+        self.n_components = n_components
+        self.metric = metric
+        self.metric_params = metric_params
+
+    def __sklearn_tags__(self):
+        tags = super().__sklearn_tags__()
+        tags.input_tags.pairwise = self.metric == "precomputed"
+        return tags
+
+    def fit(self, X, y=None):
+        """
+        Compute the embedding positions.
+
+        Parameters
+        ----------
+        X : array-like of shape (n_samples, n_features) or \
+                (n_samples, n_samples)
+            Input data. If ``metric=='precomputed'``, the input should
+            be the dissimilarity matrix.
+
+        y : Ignored
+            Not used, present for API consistency by convention.
+
+        Returns
+        -------
+        self : object
+            Fitted estimator.
+        """
+        self.fit_transform(X)
+        return self
+
+    @_fit_context(prefer_skip_nested_validation=True)
+    def fit_transform(self, X, y=None):
+        """
+        Compute and return the embedding positions.
+
+        Parameters
+        ----------
+        X : array-like of shape (n_samples, n_features) or \
+                (n_samples, n_samples)
+            Input data. If ``metric=='precomputed'``, the input should
+            be the dissimilarity matrix.
+
+        y : Ignored
+            Not used, present for API consistency by convention.
+
+        Returns
+        -------
+        X_new : ndarray of shape (n_samples, n_components)
+            The embedding coordinates.
+        """
+
+        X = validate_data(self, X)
+
+        if self.metric == "precomputed":
+            self.dissimilarity_matrix_ = X
+            self.dissimilarity_matrix_ = check_symmetric(
+                self.dissimilarity_matrix_, raise_exception=True
+            )
+        else:
+            self.dissimilarity_matrix_ = pairwise_distances(
+                X,
+                metric=self.metric,
+                **(self.metric_params if self.metric_params is not None else {}),
+            )
+
+        # Double centering
+        B = self.dissimilarity_matrix_**2
+        B = B.astype(np.float64)
+        B -= np.mean(B, axis=0)
+        B -= np.mean(B, axis=1, keepdims=True)
+        B *= -0.5
+
+        # Eigendecomposition
+        w, U = linalg.eigh(B)
+
+        # Reversing the order of the eigenvalues/eigenvectors to put
+        # the eigenvalues in decreasing order
+        w = w[::-1][: self.n_components]
+        U = U[:, ::-1][:, : self.n_components]
+
+        # Set the signs of eigenvectors to enforce deterministic output
+        U, _ = svd_flip(U, None)
+
+        self.embedding_ = np.sqrt(w) * U
+        self.eigenvalues_ = w
+
+        return self.embedding_
diff --git a/sklearn/manifold/tests/test_classical_mds.py b/sklearn/manifold/tests/test_classical_mds.py
@@ -0,0 +1,68 @@
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+from sklearn.datasets import load_iris
+from sklearn.decomposition import PCA
+from sklearn.manifold import ClassicalMDS
+from sklearn.metrics import euclidean_distances
+
+
+def test_classical_mds_equivalent_to_pca():
+    X, _ = load_iris(return_X_y=True)
+
+    cmds = ClassicalMDS(n_components=2, metric="euclidean")
+    pca = PCA(n_components=2)
+
+    Z1 = cmds.fit_transform(X)
+    Z2 = pca.fit_transform(X)
+
+    # Swap the signs if necessary
+    for comp in range(2):
+        if Z1[0, comp] < 0 and Z2[0, comp] > 0:
+            Z2[:, comp] *= -1
+
+    assert_allclose(Z1, Z2)
+
+    assert_allclose(np.sqrt(cmds.eigenvalues_), pca.singular_values_)
+
+
+def test_classical_mds_equivalent_on_data_and_distances():
+    X, _ = load_iris(return_X_y=True)
+
+    cmds = ClassicalMDS(n_components=2, metric="euclidean")
+    Z1 = cmds.fit_transform(X)
+
+    cmds = ClassicalMDS(n_components=2, metric="precomputed")
+    Z2 = cmds.fit_transform(euclidean_distances(X))
+
+    assert_allclose(Z1, Z2)
+
+
+def test_classical_mds_wrong_inputs():
+    # Non-symmetric input
+    dissim = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
+    with pytest.raises(ValueError, match="Array must be symmetric"):
+        ClassicalMDS(metric="precomputed").fit(dissim)
+
+    # Non-square input
+    dissim = np.array([[0, 1, 2], [3, 4, 5]])
+    with pytest.raises(ValueError, match="array must be 2-dimensional and square"):
+        ClassicalMDS(metric="precomputed").fit(dissim)
+
+
+def test_classical_mds_metric_params():
+    X, _ = load_iris(return_X_y=True)
+
+    cmds = ClassicalMDS(n_components=2, metric="euclidean")
+    Z1 = cmds.fit_transform(X)
+
+    cmds = ClassicalMDS(n_components=2, metric="minkowski", metric_params={"p": 2})
+    Z2 = cmds.fit_transform(X)
+
+    assert_allclose(Z1, Z2)
+
+    cmds = ClassicalMDS(n_components=2, metric="minkowski", metric_params={"p": 1})
+    Z3 = cmds.fit_transform(X)
+
+    assert not np.allclose(Z1, Z3)