Issue 1:
Right now, the content of classes_ and the shape of the output of predict_proba and decision_function are inconsistent between estimators.

If an estimator supports multioutput multiclass, classes_ will be [[0, 1]] * n_outputs and predict_proba will return an object array of length n_outputs where each is (n_samples, 2).

If the estimator supports multi-label classification, but not multioutput multiclass, classes_ will be np.arange(n_outputs) and predict_proba will be (n_samples, n_outputs).

That leads to issues such as #8773, but also runs counter to the sklearn promise of a consistent API.

Given that even basic API like the score method is broken on multioutput multiclass (see #9414), I'm not sure it's worth keeping this maintenance burden.

Alternatively we could detect in the multioutput multiclass classifiers that a problem is actually multilabel and use the appropriate shapes for that.

Either way, it's an incompatible change that's tricky to do with a deprecation cycle :-/

Issue 2:
In #19859, I am heading toward addressing some issue already reported long time ago there: #2451

We have a range of classifier supporting multilabel classification. To give a bit of feedback, I recall what a multilabel indicator matrix look like:

In [1]: from sklearn.datasets import make_multilabel_classification

In [2]: X, y = make_multilabel_classification(n_samples=10)

In [3]: y
Out[3]:
array([[0, 0, 0, 0, 0],
[0, 0, 0, 0, 1],
[0, 1, 0, 0, 1],
[0, 1, 0, 0, 1],
[1, 1, 0, 0, 0],
[0, 1, 1, 0, 0],
[0, 1, 0, 0, 0],
[0, 0, 1, 0, 0],
[0, 1, 0, 0, 0],
[0, 1, 1, 1, 1]])

Each column contain 0/1 entries and you can have several 1 for each row. This type of classification is a special case of multiclass-multioutput where each column will not be limited to 0/1 entries but to multiclass labels.

We have some issue with consistency in scikit-learn when calling the method predict_proba with classifiers handling these type of data.

Current situation in scikit-learn

We have a couple of classifiers that handle multilabel targets: tree-based models, neighbors, MLP, and ridge. In addition, the tree-based models and neighbors classifiers handle multiclass-multioutput.

In a multiclass-multioutput problem, classifier.predict_proba return a list of n_outputs containing arrays of shape (n_samples, n_classes). Since multilabel classification is a specific case of multiclass-multioutput, classifier.predict_proba will return a list of n_outputs containing arrays of shape (n_samples, 2). This should be an historical choice.

MLP and ridge (currently broken - see #19869) do not support multiclass-multioutput. There internal design make that it is easy to output an array of shape (n_samples, n_outputs) containing the probability of the positive class only when calling classifier.predict_proba.

Therefore we have two prediction format for the classifier. This RFC intends at taking a decision regarding this discrepancy.

Potential solutions

1. Keep the outputs as is

We can decide to keep the output as we currently do. One can argue that by supporting multiclass-multioutput it is natural to keep the extensive format list of arrays while being efficient with arrays only for classifiers supporting multilabel-indicator targets.

It would require to make the documentation to be explicit.

2. Make output consistent

2.1 By adding an estimator parameter

We could add an estimator parameter allowing to switch from one representation to another.

2.2 By adding a parameter to the predict_proba method

We could add a parameter to the predict_proba method to switch the output of the representation. It might still imply a classifier to store which type of target where fitted to distinguish between multiclass-multioutput and multilabel-indicator at predict.

2.3 By forcing all estimators to a specific format

In addition, we could enforce a specific format such that by default all classifiers return the same output format. In this case, it requires some deprecation cycles.