Repository: scikit-learn/scikit-learn
Pull request: FIX add support for multilabel classification in RidgeClassifier*
closes #19858 
closes #19882
build on the top of #19859

This PR fixes the `predict` method in `RidgeClassifierCV` and enables support for `RidgeClassifier` when the target is a multilabel indicator matrix.

The common tests are using #19859. The change in this PR should be limited to `ridge.py` and `test_ridge.py` once #19859 is merged.

Problem statement:
RidgeClassifierCV is not supporting properly multilabel-indicator
Our documentation state that `RidgeClassifierCV` is supporting `multilabel-indicator`.
However, a call to `predict`, will output an array of (n_samples,) containing `n_classes` instead of a matrix of `(n_samples, n_classes). The method `score` will fail in this case.

It looks like we should use a `LabelBinarizer` before to predict (I did not look at code thought).

RidgeClassifier should support multilabel-indicator
`RidgeClassifierCV` when using the GCV mode supports `multilabel-indicator`. `RidgeClassifier` should indeed be able to support this case. It would make thing consistent.

Hints from issue discussion:
This is only the case when `cv=None`
Hello! I can try fixing this if nobody else wants to?
You will also need to investigate a couple of things because I find it fishy that we are supposed to provide support for `multilabel-indicator` for `RegressionClassifierCV` with `cv=None` and not provide the same support for `RegressionClassifier`
My bad `cv=None` is the using the `gcv_mode`.
@hxyue1 You would like to take this issue intro consideration. It will be indeed easier to deal with due to my ongoing work in #19859.

take
> @hxyue1 You would like to take this issue intro consideration. It will be indeed easier to deal with due to my ongoing work in #19859.

Sorry, I'm not following...?
Ups sorry, I did a typo. I would like to take this issue as part of solving the consistency with `multilabel-indicator` common test.

Commit message: update whats new

Patch to review:
diff --git a/doc/whats_new/v1.0.rst b/doc/whats_new/v1.0.rst
index ba3f6d6d1110d..b011e203ce357 100644
--- a/doc/whats_new/v1.0.rst
+++ b/doc/whats_new/v1.0.rst
@@ -229,6 +229,11 @@ Changelog
   is now faster. This is especially noticeable on large sparse input.
   :pr:`19734` by :user:`Fred Robinson <frrad>`.
 
+- |Fix| Fix a bug in :class:`linear_model.RidgeClassifierCV` where the method
+  `predict` was performing an `argmax` on the scores obtain from
+  `decision_function` instead of returning the multilabel indicator matrix.
+  :pr:`19869` by :user:`Guillaume Lemaitre <glemaitre>`.
+
 :mod:`sklearn.manifold`
 .......................
 
@@ -276,7 +281,7 @@ Changelog
   :pr:`18649` by `Leandro Hermida <hermidalc>` and
   `Rodion Martynov <marrodion>`.
 
-- |Fix| The `fit` method of the successive halving parameter search 
+- |Fix| The `fit` method of the successive halving parameter search
   (:class:`model_selection.HalvingGridSearchCV`, and
   :class:`model_selection.HalvingRandomSearchCV`) now correctly handles the
   `groups` parameter. :pr:`19847` by :user:`Xiaoyu Chai <xiaoyuchai>`.
diff --git a/sklearn/base.py b/sklearn/base.py
index ec264b0cf5edc..be4e9a540a79b 100644
--- a/sklearn/base.py
+++ b/sklearn/base.py
@@ -767,6 +767,13 @@ class MetaEstimatorMixin:
     """Mixin class for all meta estimators in scikit-learn."""
 
 
+class MultiLabelMixin:
+    """Mixin to mark estimators that support multilabel-indicator target."""
+
+    def _more_tags(self):
+        return {"multilabel": True}
+
+
 class MultiOutputMixin:
     """Mixin to mark estimators that support multioutput."""
     def _more_tags(self):
diff --git a/sklearn/ensemble/_forest.py b/sklearn/ensemble/_forest.py
index 140c1c93e8eef..e8e4b10263ba7 100644
--- a/sklearn/ensemble/_forest.py
+++ b/sklearn/ensemble/_forest.py
@@ -51,7 +51,12 @@ class calls the ``fit`` method of each sub-estimator on random samples
 from joblib import Parallel
 
 from ..base import is_classifier
-from ..base import ClassifierMixin, RegressorMixin, MultiOutputMixin
+from ..base import (
+    ClassifierMixin,
+    MultiLabelMixin,
+    MultiOutputMixin,
+    RegressorMixin,
+)
 from ..metrics import accuracy_score, r2_score
 from ..preprocessing import OneHotEncoder
 from ..tree import (DecisionTreeClassifier, DecisionTreeRegressor,
@@ -572,7 +577,9 @@ def _accumulate_prediction(predict, X, out, lock):
                 out[i] += prediction[i]
 
 
-class ForestClassifier(ClassifierMixin, BaseForest, metaclass=ABCMeta):
+class ForestClassifier(
+    ClassifierMixin, MultiLabelMixin, BaseForest, metaclass=ABCMeta,
+):
     """
     Base class for forest of trees-based classifiers.
 
diff --git a/sklearn/linear_model/_base.py b/sklearn/linear_model/_base.py
index c80c2db622921..653e3810fea78 100644
--- a/sklearn/linear_model/_base.py
+++ b/sklearn/linear_model/_base.py
@@ -387,15 +387,15 @@ def decision_function(self, X):
 
         Parameters
         ----------
-        X : array-like or sparse matrix, shape (n_samples, n_features)
-            Samples.
+        X : {array-like, sparse matrix} of shape (n_samples, n_features)
+            The data matrix for which we want to get the confidence scores.
 
         Returns
         -------
-        array, shape=(n_samples,) if n_classes == 2 else (n_samples, n_classes)
-            Confidence scores per (sample, class) combination. In the binary
-            case, confidence score for self.classes_[1] where >0 means this
-            class would be predicted.
+        scores : ndarray of shape (n_samples,) or (n_samples, n_classes)
+            Confidence scores per `(n_sample, n_classes)` combination. In the
+            binary case, confidence score for `self.classes_[1]` where >0 means
+            this class would be predicted.
         """
         check_is_fitted(self)
 
@@ -410,13 +410,13 @@ def predict(self, X):
 
         Parameters
         ----------
-        X : array-like or sparse matrix, shape (n_samples, n_features)
-            Samples.
+        X : {array-like, sparse matrix} of shape (n_samples, n_features)
+            The data matrix for which we want to get the predictions.
 
         Returns
         -------
-        C : array, shape [n_samples]
-            Predicted class label per sample.
+        y_pred : ndarray of shape (n_samples,)
+            Vector containing the class labels for each sample.
         """
         scores = self.decision_function(X)
         if len(scores.shape) == 1:
diff --git a/sklearn/linear_model/_ridge.py b/sklearn/linear_model/_ridge.py
index 343bc6a170c9b..65cd8f7352826 100644
--- a/sklearn/linear_model/_ridge.py
+++ b/sklearn/linear_model/_ridge.py
@@ -20,7 +20,12 @@
 from ._base import LinearClassifierMixin, LinearModel
 from ._base import _deprecate_normalize, _rescale_data
 from ._sag import sag_solver
-from ..base import RegressorMixin, MultiOutputMixin, is_classifier
+from ..base import (
+    MultiLabelMixin,
+    MultiOutputMixin,
+    RegressorMixin,
+    is_classifier,
+)
 from ..utils.extmath import safe_sparse_dot
 from ..utils.extmath import row_norms
 from ..utils import check_array
@@ -1817,7 +1822,7 @@ class RidgeCV(MultiOutputMixin, RegressorMixin, _BaseRidgeCV):
     """
 
 
-class RidgeClassifierCV(LinearClassifierMixin, _BaseRidgeCV):
+class RidgeClassifierCV(LinearClassifierMixin, MultiLabelMixin, _BaseRidgeCV):
     """Ridge classifier with built-in cross-validation.
 
     See glossary entry for :term:`cross-validation estimator`.
@@ -1988,10 +1993,40 @@ def fit(self, X, y, sample_weight=None):
     def classes_(self):
         return self._label_binarizer.classes_
 
+    def predict(self, X):
+        """Predict class labels for samples in `X`.
+
+        Parameters
+        ----------
+        X : {array-like, spare matrix} of shape (n_samples, n_features)
+            The data matrix for which we want to predict the targets.
+
+        Returns
+        -------
+        y_pred : ndarray of shape (n_samples,) or (n_samples, n_outputs)
+            Vector or matrix containing the predictions. In binary and
+            multiclass problems, this is a vector containing `n_samples`. In
+            a multilabel problem, it returns a matrix of shape
+            `(n_samples, n_outputs)`.
+        """
+        if self._label_binarizer.y_type_.startswith("multilabel"):
+            # Threshold such that the negative label is -1 and positive label
+            # is 1 to use the inverse transform of the label binarizer fitted
+            # during fit.
+            scores = 2 * (self.decision_function(X) > 0) - 1
+            return self._label_binarizer.inverse_transform(scores)
+        return super().predict(X)
+
     def _more_tags(self):
         return {
             '_xfail_checks': {
                 'check_sample_weights_invariance':
                 'zero sample_weight is not equivalent to removing samples',
-            }
+                # FIXME: see
+                # https://github.com/scikit-learn/scikit-learn/issues/19858
+                # to track progress to resolve this issue
+                # 'check_classifiers_multilabel_format_output':
+                # 'RidgeClassifierCV.predict output an array of shape (25,) '
+                # 'instead of (25, 5)',
+            },
         }
diff --git a/sklearn/linear_model/tests/test_ridge.py b/sklearn/linear_model/tests/test_ridge.py
index b812788239b14..83a1445c9434d 100644
--- a/sklearn/linear_model/tests/test_ridge.py
+++ b/sklearn/linear_model/tests/test_ridge.py
@@ -1412,3 +1412,10 @@ def test_ridge_sag_with_X_fortran():
     X = X[::2, :]
     y = y[::2]
     Ridge(solver='sag').fit(X, y)
+
+
+def test_ridge_xxx():
+    # X, y = datasets.load_iris(return_X_y=True)
+    X, y = datasets.make_multilabel_classification()
+    clf = RidgeClassifierCV().fit(X, y)
+    print(clf.predict(X))
diff --git a/sklearn/neighbors/_classification.py b/sklearn/neighbors/_classification.py
index 71b869977f6aa..bf76ed5bd22b9 100644
--- a/sklearn/neighbors/_classification.py
+++ b/sklearn/neighbors/_classification.py
@@ -16,12 +16,13 @@
 import warnings
 from ._base import _check_weights, _get_weights
 from ._base import NeighborsBase, KNeighborsMixin, RadiusNeighborsMixin
-from ..base import ClassifierMixin
+from ..base import ClassifierMixin, MultiLabelMixin
 from ..utils.validation import _deprecate_positional_args
 
 
 class KNeighborsClassifier(KNeighborsMixin,
                            ClassifierMixin,
+                           MultiLabelMixin,
                            NeighborsBase):
     """Classifier implementing the k-nearest neighbors vote.
 
@@ -276,6 +277,7 @@ def predict_proba(self, X):
 
 class RadiusNeighborsClassifier(RadiusNeighborsMixin,
                                 ClassifierMixin,
+                                MultiLabelMixin,
                                 NeighborsBase):
     """Classifier implementing a vote among neighbors within a given radius
 
diff --git a/sklearn/neural_network/_multilayer_perceptron.py b/sklearn/neural_network/_multilayer_perceptron.py
index 04822360791e7..bea82a26eccd6 100644
--- a/sklearn/neural_network/_multilayer_perceptron.py
+++ b/sklearn/neural_network/_multilayer_perceptron.py
@@ -13,7 +13,12 @@
 
 import scipy.optimize
 
-from ..base import BaseEstimator, ClassifierMixin, RegressorMixin
+from ..base import (
+    BaseEstimator,
+    ClassifierMixin,
+    MultiLabelMixin,
+    RegressorMixin,
+)
 from ..base import is_classifier
 from ._base import ACTIVATIONS, DERIVATIVES, LOSS_FUNCTIONS
 from ._stochastic_optimizers import SGDOptimizer, AdamOptimizer
@@ -698,7 +703,9 @@ def _partial_fit(self, X, y):
         return self._fit(X, y, incremental=True)
 
 
-class MLPClassifier(ClassifierMixin, BaseMultilayerPerceptron):
+class MLPClassifier(
+    ClassifierMixin, MultiLabelMixin, BaseMultilayerPerceptron
+):
     """Multi-layer Perceptron classifier.
 
     This model optimizes the log-loss function using LBFGS or stochastic
diff --git a/sklearn/tree/_classes.py b/sklearn/tree/_classes.py
index 420292881f7db..5d0fd67cdcf0f 100644
--- a/sklearn/tree/_classes.py
+++ b/sklearn/tree/_classes.py
@@ -29,6 +29,7 @@
 from ..base import clone
 from ..base import RegressorMixin
 from ..base import is_classifier
+from ..base import MultiLabelMixin
 from ..base import MultiOutputMixin
 from ..utils import Bunch
 from ..utils import check_random_state
@@ -619,7 +620,9 @@ def feature_importances_(self):
 # Public estimators
 # =============================================================================
 
-class DecisionTreeClassifier(ClassifierMixin, BaseDecisionTree):
+class DecisionTreeClassifier(
+    ClassifierMixin, MultiLabelMixin, BaseDecisionTree
+):
     """A decision tree classifier.
 
     Read more in the :ref:`User Guide <tree>`.
diff --git a/sklearn/utils/estimator_checks.py b/sklearn/utils/estimator_checks.py
index 71f5b3b42de42..d0d30a2647660 100644
--- a/sklearn/utils/estimator_checks.py
+++ b/sklearn/utils/estimator_checks.py
@@ -19,6 +19,7 @@
 from ._testing import assert_array_almost_equal
 from ._testing import assert_allclose
 from ._testing import assert_allclose_dense_sparse
+from ._testing import assert_array_less
 from ._testing import set_random_state
 from ._testing import SkipTest
 from ._testing import ignore_warnings
@@ -120,6 +121,7 @@ def _yield_checks(estimator):
 
     yield check_estimator_get_tags_default_keys
 
+
 def _yield_classifier_checks(classifier):
     tags = _safe_tags(classifier)
 
@@ -139,6 +141,7 @@ def _yield_classifier_checks(classifier):
     yield check_classifiers_regression_target
     if tags["multilabel"]:
         yield check_classifiers_multilabel_representation_invariance
+        yield check_classifiers_multilabel_format_output
     if not tags["no_validation"]:
         yield check_supervised_y_no_nan
         if not tags['multioutput_only']:
@@ -637,7 +640,7 @@ def _set_checking_parameters(estimator):
         estimator.set_params(strategy='stratified')
 
     # Speed-up by reducing the number of CV or splits for CV estimators
-    loo_cv = ['RidgeCV']
+    loo_cv = ['RidgeCV', 'RidgeClassifierCV']
     if name not in loo_cv and hasattr(estimator, 'cv'):
         estimator.set_params(cv=3)
     if hasattr(estimator, 'n_splits'):
@@ -2084,15 +2087,15 @@ def check_outliers_train(name, estimator_orig, readonly_memmap=True):
                 estimator.fit(X)
 
 
-@ignore_warnings(category=(FutureWarning))
+@ignore_warnings(category=FutureWarning)
 def check_classifiers_multilabel_representation_invariance(
     name, classifier_orig
 ):
-
-    X, y = make_multilabel_classification(n_samples=100, n_features=20,
+    X, y = make_multilabel_classification(n_samples=100, n_features=2,
                                           n_classes=5, n_labels=3,
                                           length=50, allow_unlabeled=True,
                                           random_state=0)
+    X = scale(X)
 
     X_train, y_train = X[:80], y[:80]
     X_test = X[80:]
@@ -2120,6 +2123,119 @@ def check_classifiers_multilabel_representation_invariance(
     assert type(y_pred) == type(y_pred_list_of_lists)
 
 
+@ignore_warnings(category=FutureWarning)
+def check_classifiers_multilabel_format_output(name, classifier_orig):
+    """Check the output of the response methods for classifiers supporting
+    multilabel-indicator targets."""
+    classifier = clone(classifier_orig)
+    set_random_state(classifier)
+
+    n_samples, test_size, n_outputs = 100, 25, 5
+    X, y = make_multilabel_classification(n_samples=n_samples, n_features=2,
+                                          n_classes=n_outputs, n_labels=3,
+                                          length=50, allow_unlabeled=True,
+                                          random_state=0)
+    X = scale(X)
+
+    X_train, X_test = X[:-test_size], X[-test_size:]
+    y_train, y_test = y[:-test_size], y[-test_size:]
+    classifier.fit(X_train, y_train)
+
+    response_method_name = ["predict", "predict_proba", "decision_function"]
+    for method_name in response_method_name:
+        response_method = getattr(classifier, method_name, None)
+        if response_method is not None:
+            y_pred = response_method(X_test)
+
+            if method_name == "predict":
+                # y_pred.shape -> y_test.shape with the same dtype
+                assert isinstance(y_pred, np.ndarray), (
+                    f"{name}.{method_name} is expected to output a NumPy "
+                    f"array. Got {type(y_pred)} instead."
+                )
+                assert y_pred.shape == y_test.shape, (
+                    f"{name}.{method_name} output a NumPy array of shape "
+                    f"{y_pred.shape} instead of {y_test.shape}."
+                )
+                assert y_pred.dtype == y_test.dtype, (
+                    f"{name}.{method_name} does not output the same dtype than"
+                    f" the targets. Got {y_pred.dtype} instead of "
+                    f"{y_test.dtype}."
+                )
+            elif method_name == "predict_proba":
+                # y_pred.shape -> 2 possibilities:
+                # - list of length n_outputs of shape (n_samples, 2);
+                # - ndarray of shape (n_samples, n_outputs).
+                # dtype should be floating
+                if isinstance(y_pred, list):
+                    assert len(y_pred) == n_outputs, (
+                        f"{name}.{method_name} is expected to output a list "
+                        f"of length n_outputs of Numpy array. Got length of "
+                        f"{len(y_pred)} instead of {n_outputs}."
+                    )
+                    for pred in y_pred:
+                        assert pred.shape == (test_size, 2), (
+                            f"{name}.{method_name} is expected to output a "
+                            f"list of NumPy array of shape (n_samples, 2). Got"
+                            f" {pred.shape} instead of {(test_size, 2)}."
+                        )
+                        assert pred.dtype.kind == "f", (
+                            f"{name}.{method_name} is expected to output a "
+                            f"list of NumPy array of floating dtype. Got "
+                            f"{pred.dtype} instead."
+                        )
+                        # check that we have the correct probabilities
+                        assert_allclose(
+                            pred.sum(axis=1),
+                            1,
+                            err_msg=(
+                                f"{name}.{method_name} is expected to provide "
+                                f"probabilities such that each array rows "
+                                f"should sum to 1."
+                            ),
+                        )
+                elif isinstance(y_pred, np.ndarray):
+                    assert y_pred.shape == (test_size, n_outputs), (
+                        f"{name}.{method_name} is expected to output a NumPy "
+                        f"array of shape (n_samples, n_outputs). Got "
+                        f"{y_pred.shape} instead of {(test_size, n_outputs)}."
+                    )
+                    assert y_pred.dtype.kind == "f", (
+                        f"{name}.{method_name} is expected to output a NumPy "
+                        f"array of floating dtype. Got {y_pred.dtype} instead."
+                    )
+                    assert_array_less(
+                        y_pred,
+                        1,
+                        err_msg=(
+                            f"{name}.{method_name} is expected to provide "
+                            f"probabilities of the positive class and should "
+                            f"therefore contain values below 1."
+                        ),
+                    )
+                else:
+                    raise ValueError(
+                        f"Unknown returned type {type(y_pred)} by "
+                        f"{name}.{method_name}. A list or a Numpy array are "
+                        f"expected."
+                    )
+            else:  # "decision_function"
+                # y_pred.shape -> y_test.shape with floating dtype
+                assert isinstance(y_pred, np.ndarray), (
+                    f"{name}.{method_name} is expected to output a NumPy "
+                    f"array. Got {type(y_pred)} instead."
+                )
+                assert y_pred.shape == (test_size, n_outputs), (
+                    f"{name}.{method_name} is expected to provide a NumPy "
+                    f"array of shape (n_samples, n_outputs). Got "
+                    f"{y_pred.shape} instead of {(test_size, n_outputs)}."
+                )
+                assert y_pred.dtype.kind == "f", (
+                    f"{name}.{method_name} is expected to output a floating "
+                    f"dtype. Got {y_pred.dtype} instead."
+                )
+
+
 @ignore_warnings(category=FutureWarning)
 def check_estimators_fit_returns_self(
     name, estimator_orig, readonly_memmap=False
diff --git a/sklearn/utils/tests/test_estimator_checks.py b/sklearn/utils/tests/test_estimator_checks.py
index 4792f50f2baef..f3b92c3678aff 100644
--- a/sklearn/utils/tests/test_estimator_checks.py
+++ b/sklearn/utils/tests/test_estimator_checks.py
@@ -10,6 +10,7 @@
 import joblib
 
 from sklearn.base import BaseEstimator, ClassifierMixin
+from sklearn.datasets import make_multilabel_classification
 from sklearn.utils import deprecated
 from sklearn.utils._testing import (
     assert_raises,
@@ -21,18 +22,6 @@
     MinimalTransformer,
     SkipTest,
 )
-from sklearn.utils.estimator_checks import check_estimator, _NotAnArray
-from sklearn.utils.estimator_checks \
-    import check_class_weight_balanced_linear_classifier
-from sklearn.utils.estimator_checks import set_random_state
-from sklearn.utils.estimator_checks import _set_checking_parameters
-from sklearn.utils.estimator_checks import check_estimators_unfitted
-from sklearn.utils.estimator_checks import check_fit_score_takes_y
-from sklearn.utils.estimator_checks import check_no_attributes_set_in_init
-from sklearn.utils.estimator_checks import check_classifier_data_not_an_array
-from sklearn.utils.estimator_checks import check_regressor_data_not_an_array
-from sklearn.utils.estimator_checks import \
-    check_estimator_get_tags_default_keys
 from sklearn.utils.validation import check_is_fitted
 from sklearn.utils.estimator_checks import check_outlier_corruption
 from sklearn.utils.fixes import np_version, parse_version
@@ -48,6 +37,21 @@
 from sklearn.utils import all_estimators
 from sklearn.exceptions import SkipTestWarning
 
+from sklearn.utils.estimator_checks import (
+    _NotAnArray,
+    _set_checking_parameters,
+    check_class_weight_balanced_linear_classifier,
+    check_classifier_data_not_an_array,
+    check_classifiers_multilabel_format_output,
+    check_estimator,
+    check_estimator_get_tags_default_keys,
+    check_estimators_unfitted,
+    check_fit_score_takes_y,
+    check_no_attributes_set_in_init,
+    check_regressor_data_not_an_array,
+    set_random_state,
+)
+
 
 class CorrectNotFittedError(ValueError):
     """Exception class to raise if estimator is used before fitting.
@@ -673,6 +677,244 @@ def test_check_estimator_get_tags_default_keys():
     )
 
 
+def test_check_classifiers_multilabel_output_format():
+    n_samples, test_size, n_outputs = 100, 25, 5
+    _, y = make_multilabel_classification(n_samples=n_samples, n_features=2,
+                                          n_classes=n_outputs, n_labels=3,
+                                          length=50, allow_unlabeled=True,
+                                          random_state=0)
+    y_test = y[-test_size:]
+
+    class BaseMultiLabelClassifierMock(ClassifierMixin, BaseEstimator):
+        def __init__(self, response_output):
+            self.response_output = response_output
+
+        def fit(self, X, y):
+            return self
+
+    class MultiLabelClassifierPredict(BaseMultiLabelClassifierMock):
+        def predict(self, X):
+            return self.response_output
+
+    class MultiLabelClassifierPredictProba(BaseMultiLabelClassifierMock):
+        def predict_proba(self, X):
+            return self.response_output
+
+    class MultiLabelClassifierDecisionFunction(BaseMultiLabelClassifierMock):
+        def decision_function(self, X):
+            return self.response_output
+
+    # 1.method predict
+    # 1.1 inconsistent array type
+    clf = MultiLabelClassifierPredict(response_output=y_test.tolist())
+    err_msg = (
+        r"MultiLabelClassifierPredict.predict is expected to output a "
+        r"NumPy array. Got <class 'list'> instead."
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+    # 1.2. inconsistent shape
+    clf = MultiLabelClassifierPredict(response_output=y_test[:, :-1])
+    err_msg = (
+        r"MultiLabelClassifierPredict.predict output a NumPy array of "
+        r"shape \(25, 4\) instead of \(25, 5\)."
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+    # 1.3 inconsistent dtype
+    clf = MultiLabelClassifierPredict(
+        response_output=y_test.astype(np.float64)
+    )
+    err_msg = (
+        r"MultiLabelClassifierPredict.predict does not output the same "
+        r"dtype than the targets."
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+
+    # 2. method predict_proba
+    # 2.1 unknown output type
+    clf = MultiLabelClassifierPredictProba(
+        response_output=sp.csr_matrix(y_test)
+    )
+    err_msg = (
+        r"Unknown returned type <class 'scipy.sparse.csr.csr_matrix'> by "
+        r"MultiLabelClassifierPredictProba.predict_proba. A list or a Numpy "
+        r"array are expected."
+    )
+    assert_raises_regex(
+        ValueError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+    # 2.2 for list output
+    # 2.2.1 iconsistent length
+    clf = MultiLabelClassifierPredictProba(response_output=y_test.tolist())
+    err_msg = (
+        r"MultiLabelClassifierPredictProba.predict_proba is expected to "
+        r"output a list of length n_outputs of Numpy array. Got length of 25 "
+        r"instead of 5."
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+    # 2.2.2 array of inconsistent shape
+    response_output = [np.ones_like(y_test) for _ in range(n_outputs)]
+    clf = MultiLabelClassifierPredictProba(response_output=response_output)
+    err_msg = (
+        r"MultiLabelClassifierPredictProba.predict_proba is expected to output"
+        r" a list of NumPy array of shape \(n_samples, 2\). Got \(25, 5\) "
+        r"instead of \(25, 2\)."
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+    # 2.2.3 array of inconsistent dtype
+    response_output = [
+        np.ones(shape=(y_test.shape[0], 2), dtype=np.int64)
+        for _ in range(n_outputs)
+    ]
+    clf = MultiLabelClassifierPredictProba(response_output=response_output)
+    err_msg = (
+        r"MultiLabelClassifierPredictProba.predict_proba is expected to output"
+        r" a list of NumPy array of floating dtype."
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+    # 2.2.4 array does not contain probability (each row should sum to 1)
+    response_output = [
+        np.ones(shape=(y_test.shape[0], 2), dtype=np.float64)
+        for _ in range(n_outputs)
+    ]
+    clf = MultiLabelClassifierPredictProba(response_output=response_output)
+    err_msg = (
+        r"MultiLabelClassifierPredictProba.predict_proba is expected to "
+        r"provide probabilities such that each array rows should sum to 1."
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+    # 2.3 for array output
+    # 2.3.1 array of inconsistent shape
+    clf = MultiLabelClassifierPredictProba(response_output=y_test[:, :-1])
+    err_msg = (
+        r"MultiLabelClassifierPredictProba.predict_proba is expected to "
+        r"output a NumPy array of shape \(n_samples, n_outputs\). Got "
+        r"\(25, 4\) instead of \(25, 5\)."
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+    # 2.3.2 array of inconsistent dtype
+    response_output = np.zeros_like(y_test, dtype=np.int64)
+    clf = MultiLabelClassifierPredictProba(response_output=response_output)
+    err_msg = (
+        r"MultiLabelClassifierPredictProba.predict_proba is expected to "
+        r"output a NumPy array of floating dtype."
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+    # 2.2.4 array does not contain probabilities
+    clf = MultiLabelClassifierPredictProba(response_output=y_test * 2.0)
+    err_msg = (
+        r"MultiLabelClassifierPredictProba.predict_proba is expected to "
+        r"provide probabilities of the positive class and should therefore "
+        r"contain values below 1."
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+
+    # 3. decision_function
+    # 3.1 inconsistent array type
+    clf = MultiLabelClassifierDecisionFunction(response_output=y_test.tolist())
+    err_msg = (
+        r"MultiLabelClassifierDecisionFunction.decision_function is expected "
+        r"to output a NumPy array. Got <class 'list'> instead."
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+    # 3.2. inconsistent shape
+    clf = MultiLabelClassifierDecisionFunction(response_output=y_test[:, :-1])
+    err_msg = (
+        r"MultiLabelClassifierDecisionFunction.decision_function is expected "
+        r"to provide a NumPy array of shape \(n_samples, n_outputs\). Got "
+        r"\(25, 4\) instead of \(25, 5\)"
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+    # 3.3 inconsistent dtype
+    clf = MultiLabelClassifierDecisionFunction(response_output=y_test)
+    err_msg = (
+        r"MultiLabelClassifierDecisionFunction.decision_function is expected "
+        r"to output a floating dtype."
+    )
+    assert_raises_regex(
+        AssertionError,
+        err_msg,
+        check_classifiers_multilabel_format_output,
+        clf.__class__.__name__,
+        clf,
+    )
+
+
 def run_tests_without_pytest():
     """Runs the tests in this file without using pytest.
     """


Code review comments:
diff_hunk: @@ -2120,6 +2123,119 @@ def check_classifiers_multilabel_representation_invariance(
     assert type(y_pred) == type(y_pred_list_of_lists)
 
 
+@ignore_warnings(category=FutureWarning)
+def check_classifiers_multilabel_format_output(name, classifier_orig):
+    """Check the output of the response methods for classifiers supporting
+    multilabel-indicator targets."""
+    classifier = clone(classifier_orig)
+    set_random_state(classifier)
+
+    n_samples, test_size, n_outputs = 100, 25, 5
+    X, y = make_multilabel_classification(n_samples=n_samples, n_features=2,
+                                          n_classes=n_outputs, n_labels=3,
+                                          length=50, allow_unlabeled=True,
+                                          random_state=0)
+    X = scale(X)
+
+    X_train, X_test = X[:-test_size], X[-test_size:]
+    y_train, y_test = y[:-test_size], y[-test_size:]
+    classifier.fit(X_train, y_train)
+
+    response_method_name = ["predict", "predict_proba", "decision_function"]
+    for method_name in response_method_name:
+        response_method = getattr(classifier, method_name, None)
+        if response_method is not None:
+            y_pred = response_method(X_test)
+
+            if method_name == "predict":
path: sklearn/utils/estimator_checks.py
text: @user1:
Since `method_name` all has its own distinct path, I think we do not need the for loop? 

```
predict_method = getattr(classifier, "predict", None)
if predict_method is not None:
    ...

predict_proba_method = ...
if predict_proba_method is not None:
    ...
```

This would also help with lower the indentations.

@author:
OK. Be aware that this comment will be addressed in #19859 :)
I did not want to write specific test for this issue knowing that it should be treated in the common tests.