Adding multioutput support for KPCovC

Author: rvasav26

src/skmatter/decomposition/_kernel_pcovc.py

@@ -1,6 +1,7 @@
 import numpy as np
 
 from sklearn import clone
+from sklearn.multioutput import MultiOutputClassifier
 from sklearn.svm import LinearSVC
 from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
 from sklearn.linear_model import (
@@ -52,6 +53,9 @@ class KernelPCovC(LinearClassifierMixin, _BaseKPCov):
 
             n_components == n_samples
 
+    n_outputs : int
+        The number of outputs when ``fit`` is performed.
+
     svd_solver : {'auto', 'full', 'arpack', 'randomized'}, default='auto'
         If auto :
             The solver is selected by a default policy based on `X.shape` and
@@ -78,13 +82,14 @@ class KernelPCovC(LinearClassifierMixin, _BaseKPCov):
         - ``sklearn.linear_model.LogisticRegressionCV()``
         - ``sklearn.svm.LinearSVC()``
         - ``sklearn.discriminant_analysis.LinearDiscriminantAnalysis()``
+        - ``sklearn.multioutput.MultiOutputClassifier()``
         - ``sklearn.linear_model.RidgeClassifier()``
         - ``sklearn.linear_model.RidgeClassifierCV()``
         - ``sklearn.linear_model.Perceptron()``
 
         If a pre-fitted classifier is provided, it is used to compute :math:`{\mathbf{Z}}`.
-        If None, ``sklearn.linear_model.LogisticRegression()``
-        is used as the classifier.
+        If None and ``n_outputs < 2``, ``sklearn.linear_model.LogisticRegression()`` is used.
+        If None and ``n_outputs == 2``, ``sklearn.multioutput.MultiOutputClassifier()`` is used.
 
     kernel : {"linear", "poly", "rbf", "sigmoid", "precomputed"} or callable, default="linear"
         Kernel.
@@ -132,6 +137,9 @@ class KernelPCovC(LinearClassifierMixin, _BaseKPCov):
 
     Attributes
     ----------
+    n_outputs : int
+        The number of outputs when ``fit`` is performed.
+
     classifier : estimator object
         The linear classifier passed for fitting. If pre-fitted, it is assummed
         to be fit on a precomputed kernel :math:`\mathbf{K}` and :math:`\mathbf{Y}`.
@@ -268,9 +276,11 @@ def fit(self, X, Y, W=None):
         self: object
             Returns the instance itself.
         """
-        X, Y = validate_data(self, X, Y, y_numeric=False)
+        X, Y = validate_data(self, X, Y, multi_output=True, y_numeric=False)
+
         check_classification_targets(Y)
         self.classes_ = np.unique(Y)
+        self.n_outputs = 1 if Y.ndim == 1 else Y.shape[1]
 
         super()._set_fit_params(X)
 
@@ -285,6 +295,7 @@ def fit(self, X, Y, W=None):
             LogisticRegressionCV,
             LinearSVC,
             LinearDiscriminantAnalysis,
+            MultiOutputClassifier,
             RidgeClassifier,
             RidgeClassifierCV,
             SGDClassifier,
@@ -300,27 +311,37 @@ def fit(self, X, Y, W=None):
                 ", or `precomputed`"
             )
 
-        if self.classifier != "precomputed":
-            if self.classifier is None:
-                classifier = LogisticRegression()
-            else:
-                classifier = self.classifier
+        multioutput = self.n_outputs != 1
+        precomputed = self.classifier == "precomputed"
 
-            # for convergence warnings
-            if hasattr(classifier, "max_iter") and (
-                classifier.max_iter is None or classifier.max_iter < 500
-            ):
-                classifier.max_iter = 500
+        if self.classifier is None or precomputed:
+            # used as the default classifier for subsequent computations
+            classifier = (
+                MultiOutputClassifier(LogisticRegression())
+                if multioutput
+                else LogisticRegression()
+            )
+        else:
+            classifier = self.classifier
 
-            # Check if classifier is fitted; if not, fit with precomputed K
-            self.z_classifier_ = check_cl_fit(classifier, K, Y)
-            W = self.z_classifier_.coef_.T
+        if hasattr(classifier, "max_iter") and (
+            classifier.max_iter is None or classifier.max_iter < 500
+        ):
+            classifier.max_iter = 500
+
+        if precomputed and W is None:
+            _ = clone(classifier).fit(K, Y)
+            if multioutput:
+                W = np.hstack([_.coef_.T for _ in _.estimators_])
+            else:
+                W = _.coef_.T
 
         else:
-            # If precomputed, use default classifier to predict Y from T
-            classifier = LogisticRegression(max_iter=500)
-            if W is None:
-                W = LogisticRegression().fit(K, Y).coef_.T
+            self.z_classifier_ = check_cl_fit(classifier, K, Y)
+            if multioutput:
+                W = np.hstack([est_.coef_.T for est_ in self.z_classifier_.estimators_])
+            else:
+                W = self.z_classifier_.coef_.T
 
         Z = K @ W
 
@@ -333,10 +354,16 @@ def fit(self, X, Y, W=None):
 
         self.classifier_ = clone(classifier).fit(K @ self.pkt_, Y)
 
-        self.ptz_ = self.classifier_.coef_.T
-        self.pkz_ = self.pkt_ @ self.ptz_
+        if multioutput:
+            self.ptz_ = np.hstack(
+                [est_.coef_.T for est_ in self.classifier_.estimators_]
+            )
+            self.pkz_ = self.pkt_ @ self.ptz_
+        else:
+            self.ptz_ = self.classifier_.coef_.T
+            self.pkz_ = self.pkt_ @ self.ptz_
 
-        if len(Y.shape) == 1 and type_of_target(Y) == "binary":
+        if not multioutput and type_of_target(Y) == "binary":
             self.pkz_ = self.pkz_.reshape(
                 K.shape[1],
             )
@@ -345,6 +372,7 @@ def fit(self, X, Y, W=None):
             )
 
         self.components_ = self.pkt_.T  # for sklearn compatibility
+
         return self
 
     def predict(self, X=None, T=None):
@@ -424,9 +452,12 @@ def decision_function(self, X=None, T=None):
 
         Returns
         -------
-        Z : numpy.ndarray, shape (n_samples,) or (n_samples, n_classes)
+        Z : numpy.ndarray, shape (n_samples,) or (n_samples, n_classes), or a list of \
+                n_outputs such arrays if n_outputs > 1
             Confidence scores. For binary classification, has shape `(n_samples,)`,
-            for multiclass classification, has shape `(n_samples, n_classes)`
+            for multiclass classification, has shape `(n_samples, n_classes)`. 
+            If n_outputs > 1, the list can contain arrays with differing shapes 
+            depending on the number of classes in each output of Y.
         """
         check_is_fitted(self, attributes=["pkz_", "ptz_"])
 
@@ -439,9 +470,21 @@ def decision_function(self, X=None, T=None):
             if self.center:
                 K = self.centerer_.transform(K)
 
-            # Or self.classifier_.decision_function(K @ self.pkt_)
-            return K @ self.pkz_ + self.classifier_.intercept_
+            if self.n_outputs == 1:
+                # Or self.classifier_.decision_function(K @ self.pkt_)
+                return K @ self.pkz_ + self.classifier_.intercept_
+            else:
+                return [
+                    est_.decision_function(K @ self.pkt_)
+                    for est_ in self.classifier_.estimators_
+                ]
 
         else:
             T = check_array(T)
-            return T @ self.ptz_ + self.classifier_.intercept_
+
+            if self.n_outputs == 1:
+                T @ self.ptz_ + self.classifier_.intercept_
+            else:
+                return [
+                    est_.decision_function(T) for est_ in self.classifier_.estimators_
+                ]

src/skmatter/decomposition/_pcovc.py

@@ -22,11 +22,11 @@
 
 
 # No inheritance from MultiOutputMixin because decision_function would fail
-# test_check_estimator.py 'check_classifier_multioutput' (line 2479 of estimator_checks.py)
-# - this is the only test for MultiOutputClassifiers, so is it OK to exclude this tag?
+# test_check_estimator.py 'check_classifier_multioutput' (line 2479 of estimator_checks.py).
+# This is the only test for multioutput classifiers, so is it OK to exclude this tag?
 
 # did a search of all classifiers that inherit from MultiOutputMixin - none of them implement
-# decision function, so I don't think we need to inherit
+# decision function
 
 
 class PCovC(LinearClassifierMixin, _BasePCov):
@@ -120,6 +120,7 @@ class PCovC(LinearClassifierMixin, _BasePCov):
         - ``sklearn.linear_model.LogisticRegressionCV()``
         - ``sklearn.svm.LinearSVC()``
         - ``sklearn.discriminant_analysis.LinearDiscriminantAnalysis()``
+        - ``sklearn.multioutput.MultiOutputClassifier()``
         - ``sklearn.linear_model.RidgeClassifier()``
         - ``sklearn.linear_model.RidgeClassifierCV()``
         - ``sklearn.linear_model.Perceptron()``
@@ -131,8 +132,8 @@ class PCovC(LinearClassifierMixin, _BasePCov):
         `sklearn.pipeline.Pipeline` with model caching.
         In such cases, the classifier will be re-fitted on the same
         training data as the composite estimator.
-        If None and ``Y.ndim < 2``, ``sklearn.linear_model.LogisticRegression()`` is used.
-        If None and ``Y.ndim == 2``, ``sklearn.multioutput.MultiOutputClassifier()`` is used.
+        If None and ``n_outputs < 2``, ``sklearn.linear_model.LogisticRegression()`` is used.
+        If None and ``n_outputs == 2``, ``sklearn.multioutput.MultiOutputClassifier()`` is used.
 
     iterated_power : int or 'auto', default='auto'
         Number of iterations for the power method computed by
@@ -164,6 +165,9 @@ class PCovC(LinearClassifierMixin, _BasePCov):
         n_components, or the lesser value of n_features and n_samples
         if n_components is None.
 
+    n_outputs : int
+        The number of outputs when ``fit`` is performed.
+
     classifier : estimator object
         The linear classifier passed for fitting.
 
@@ -263,16 +267,14 @@ def fit(self, X, Y, W=None):
 
         Y : numpy.ndarray, shape (n_samples,) or (n_samples, n_outputs)
             Training data, where n_samples is the number of samples and
-            n_outputs is the number of outputs. If ``self.classifier`` is an instance
-            of ``sklearn.multioutput.MultiOutputClassifier()``, Y can be of shape
-            (n_samples, n_outputs).
+            n_outputs is the number of outputs.
 
         W : numpy.ndarray, shape (n_features, n_classes)
             Classification weights, optional when classifier is ``precomputed``. If
             not passed, it is assumed that the weights will be taken from a
             linear classifier fit between :math:`\mathbf{X}` and :math:`\mathbf{Y}`.
-            In the case of a multioutput classifier ``classifier``,
-           `` W = np.hstack([est_.coef_.T for est_ in classifier.estimators_])``.
+            In the multioutput case,
+            `` W = np.hstack([est_.coef_.T for est_ in classifier.estimators_])``.
         """
         X, Y = validate_data(self, X, Y, multi_output=True, y_numeric=False)
 
@@ -303,49 +305,31 @@ def fit(self, X, Y, W=None):
                 ", or `precomputed`"
             )
 
-        if self.n_outputs == 1 and isinstance(self.classifier, MultiOutputClassifier):
-            raise ValueError(
-                "Classifier cannot be an instance of `MultiOutputClassifier` when Y is 1D"
-            )
+        multioutput = self.n_outputs != 1
+        precomputed = self.classifier == "precomputed"
 
-        if (
-            self.n_outputs != 1
-            and self.classifier not in ["precomputed", None]
-            and not (
-                isinstance(self.classifier, MultiOutputClassifier)
-                or self.classifier == "precomputed"
-            )
-        ):
-            raise ValueError(
-                "Classifier must be an instance of `MultiOutputClassifier` when Y is 2D"
+        if self.classifier is None or precomputed:
+            # used as the default classifier for subsequent computations
+            classifier = (
+                MultiOutputClassifier(LogisticRegression())
+                if multioutput
+                else LogisticRegression()
             )
+        else:
+            classifier = self.classifier
 
-        if self.n_outputs == 1:
-            if self.classifier != "precomputed":
-                classifier = self.classifier or LogisticRegression()
-                self.z_classifier_ = check_cl_fit(classifier, X, Y)
-                W = self.z_classifier_.coef_.T
-
+        if precomputed and W is None:
+            _ = clone(classifier).fit(X, Y)
+            if multioutput:
+                W = np.hstack([_.coef_.T for _ in _.estimators_])
             else:
-                # to be used later on as the classifier fit between T and Y
-                classifier = LogisticRegression()
-                if W is None:
-                    W = clone(classifier).fit(X, Y).coef_.T
-
+                W = _.coef_.T
         else:
-            if self.classifier != "precomputed":
-                classifier = self.classifier or MultiOutputClassifier(
-                    estimator=LogisticRegression()
-                )
-                self.z_classifier_ = check_cl_fit(classifier, X, Y)
+            self.z_classifier_ = check_cl_fit(classifier, X, Y)
+            if multioutput:
                 W = np.hstack([est_.coef_.T for est_ in self.z_classifier_.estimators_])
-
             else:
-                # to be used later on as the classifier fit between T and Y
-                classifier = MultiOutputClassifier(estimator=LogisticRegression())
-                if W is None:
-                    _ = clone(classifier).fit(X, Y)
-                    W = np.hstack([_.coef_.T for _ in _.estimators_])
+                W = self.z_classifier_.coef_.T
 
         Z = X @ W
 
@@ -358,21 +342,21 @@ def fit(self, X, Y, W=None):
         # classifier and steal weights to get pxz and ptz
         self.classifier_ = clone(classifier).fit(X @ self.pxt_, Y)
 
-        if self.n_outputs == 1:
-            self.ptz_ = self.classifier_.coef_.T
-            # print(self.ptz_.shape)
-            self.pxz_ = self.pxt_ @ self.ptz_
-        else:
+        if multioutput:
             self.ptz_ = np.hstack(
                 [est_.coef_.T for est_ in self.classifier_.estimators_]
             )
             # print(f"pxt {self.pxt_.shape}")
             # print(f"ptz {self.ptz_.shape}")
             self.pxz_ = self.pxt_ @ self.ptz_
             # print(f"pxz {self.pxz_.shape}")
+        else:
+            self.ptz_ = self.classifier_.coef_.T
+            # print(self.ptz_.shape)
+            self.pxz_ = self.pxt_ @ self.ptz_
 
         # print(self.ptz_.shape)
-        if len(Y.shape) == 1 and type_of_target(Y) == "binary":
+        if not multioutput and type_of_target(Y) == "binary":
             self.pxz_ = self.pxz_.reshape(
                 X.shape[1],
             )
@@ -472,9 +456,9 @@ def decision_function(self, X=None, T=None):
         Z : numpy.ndarray, shape (n_samples,) or (n_samples, n_classes), or a list of \
                 n_outputs such arrays if n_outputs > 1
             Confidence scores. For binary classification, has shape `(n_samples,)`,
-            for multiclass classification, has shape `(n_samples, n_classes)`. If n_outputs > 1,
-            the list can contain arrays with differing shapes depending on the
-            number of classes in each output of Y.
+            for multiclass classification, has shape `(n_samples, n_classes)`. 
+            If n_outputs > 1, the list can contain arrays with differing shapes 
+            depending on the number of classes in each output of Y.
         """
         check_is_fitted(self, attributes=["pxz_", "ptz_"])
 
@@ -529,36 +513,3 @@ def transform(self, X=None):
             and n_features is the number of features.
         """
         return super().transform(X)
-
-    # def score(self, X, Y, sample_weight=None):
-    #     """Return the accuracy on the given test data and labels. Contains support
-    #     for multiclass-multioutput data.
-
-    #     Parameters
-    #     ----------
-    #     X : array-like of shape (n_samples, n_features)
-    #         Test samples.
-
-    #     Y : array-like of shape (n_samples,) or (n_samples, n_outputs)
-    #         True labels for `X`.
-
-    #     sample_weight : array-like of shape (n_samples,), default=None
-    #         Sample weights. Can only be used if the PCovC instance
-    #         has been trained on single-target data.
-
-    #     Returns
-    #     -------
-    #     score : float
-    #         Accuracy scores. If the PCovC instance was trained on a 1D Y,
-    #         this will call the ``score()`` function defined by
-    #         ``sklearn.base.ClassifierMixin``. If trained on a 2D Y, this will
-    #         call the ``score()`` function defined by
-    #         ``sklearn.multioutput.MultiOutputClassifier``.
-    #     """
-    #     X, Y = validate_data(self, X, Y, reset=False)
-
-    #     if isinstance(self.classifier_, MultiOutputClassifier):
-    #         # LinearClassifierMixin.score fails with multioutput-multiclass Y
-    #         return self.classifier_.score(X @ self.pxt_, Y)
-    #     else:
-    #         return self.classifier_.score(X @ self.pxt_, Y, sample_weight=sample_weight)