[RFC] Varying the number of outputs considered for splitting in Multi Output Decision Trees

### Describe the workflow you want to enable

One strength of RFRs is that they are incredibly robust and therefore provide a strong baseline for many tasks without needing to consider normalization or scaling of either the inputs or outputs. In the case of multi-output RFRs this robustness towards the output space goes away due to the summing of impurities across different output dimensions which entails the need to standardize the output labels to ensure that undue attention isn't given to particular outputs. Currently the documentation doesn't readily inform the user of this artifact. In the spirit of the Random Forest one solution to avoiding this problem would be to randomly sample which output(s) to consider for the determination of the split. If the number of outputs was set to 1 then we would end up with a case where the normalization of the output space once again doesn't matter.

### Describe your proposed solution

The introduction of a new kwarg `max_outputs` (by analogy to `max_features`) could allow users to control how many outputs were considered when selecting the optimal split. If set to `1.0` all outputs would be used as currently, if set to `1` then a single output would be used as described above. This seems like a relevant and natural hyper-parameter for the multi-output RFR.  

### Describe alternatives you've considered, if relevant

_No response_

### Additional context

I have not been able to find literature that explores the above slight adjustment to the current algorithm. This is a RFC to see if the team would accept such a PR in principle without the quoted 200+ citation requirement if sufficient empirical evidence was provided.

COMMENT:
I was trying to reproduce the code, was this what you had in mind? I was trying to artificially increase the MSE of y1, but I could not consistently raise it. 

@CompRhys , was this what you had in mind for reproducible code, where we try to get a significantly different value for mse_output1? Should I try more features? Thank you. 

```
import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

# Set random seed for reproducibility
np.random.seed(42)

# Generate a synthetic dataset with two features and two output dimensions
# The second output dimension has a larger scale
X = np.random.rand(100, 2)
y1 = 2 * X[:, 0] + 3 * X[:, 1] + np.random.normal(scale=0.1, size=100)
y2 = 20 * X[:, 0] + 30 * X[:, 1] + np.random.normal(scale=1, size=100)
y3 = 20 * X[:, 0] + 30 * X[:, 1] + np.random.normal(scale=1, size=100)
y4 = 20 * X[:, 0] + 30 * X[:, 1] + np.random.normal(scale=1, size=100)
y5 = 20 * X[:, 0] + 30 * X[:, 1] + np.random.normal(scale=1, size=100)
y6 = 20 * X[:, 0] + 30 * X[:, 1] + np.random.normal(scale=1, size=100)

y = np.column_stack((y1, y2,y3,y4,y5,y6))

# Split the dataset into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Create a multi-output RandomForestRegressor
rf_regressor = RandomForestRegressor(n_estimators=100, random_state=42)

# Fit the model
rf_regressor.fit(X_train, y_train)

# Make predictions
y_pred = rf_regressor.predict(X_test)

# Evaluate the model performance on each output dimension
mse_output1 = mean_squared_error(y_test[:, 0], y_pred[:, 0])
mse_output2 = mean_squared_error(y_test[:, 1], y_pred[:, 1])
mse_output3 = mean_squared_error(y_test[:, 1], y_pred[:, 1])
print("Mean Squared Error (Output 1):", mse_output1)
```
Also, while I could not find anything related specifically towards the optimization increases related to selectively choosing MSE  on each partition, I did find this paper that mentions some optimizations that could be implemented towards Random Forest. 

https://arxiv.org/pdf/2309.01460.pdf

Labels: New Feature, help wanted