[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:
Here's a slightly more involved setup to play with. I think your issue might have been that your functions are all strongly rank correlated and so the splits locations won't vary as you vary which output dominates the impurity.

```python
# %%
import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from scipy.stats import spearmanr

from sklearn.datasets import make_friedman1
from sklearn.utils.validation import check_random_state

from sklearn.compose import TransformedTargetRegressor
from sklearn.preprocessing import MinMaxScaler

# %%
# Generate a synthetic dataset with two features and two output dimensions
# The second output dimension has a larger scale
PRNG_SEED = 42
n_estimators = 100
n_samples = 100
base_noise = 0.1
f2_scale = 1
f2_noise_scale = 1
f3_scale = 10
f3_noise_scale = 10

# Generate the dataset
X, y_friedman1 = make_friedman1(n_samples=n_samples, noise=base_noise, random_state=PRNG_SEED)

generator = check_random_state(PRNG_SEED)

y_friedman2 = f2_scale * ((
        X[:, 0] ** 2 + (X[:, 1] * X[:, 2] - 1 / (X[:, 1] * X[:, 3])) ** 2
    ) ** 0.5 + f2_noise_scale * generator.standard_normal(size=(n_samples)))

y_friedman3 = f3_scale * (np.arctan(
        (X[:, 1] * X[:, 2] - 1 / (X[:, 1] * X[:, 3])) / X[:, 0]
    ) + f3_noise_scale * generator.standard_normal(size=(n_samples)))

# Check the correlation between the output dimensions
print(f"Correlation between output 1 and output 2: {spearmanr(y_friedman1, y_friedman2)}")
print(f"Correlation between output 1 and output 3: {spearmanr(y_friedman1, y_friedman3)}")
print(f"Correlation between output 2 and output 3: {spearmanr(y_friedman2, y_friedman3)}")

y = np.column_stack((y_friedman1, y_friedman2, y_friedman3))

# 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)

# %% 
rf_regressor = RandomForestRegressor(n_estimators=n_estimators, random_state=PRNG_SEED)

# %%
# Create a single-output RandomForestRegressors
# Fit the model for Output 0
rf_regressor.fit(X_train, y_train[:, 0])
y_pred = rf_regressor.predict(X_test)
mse_output1 = mean_squared_error(y_test[:, 0], y_pred)
print("Mean Squared Error (Output 0):", mse_output1)

# Fit the model for Output 1
rf_regressor.fit(X_train, y_train[:, 1])
y_pred = rf_regressor.predict(X_test)
mse_output2 = mean_squared_error(y_test[:, 1], y_pred)
print("Mean Squared Error (Output 1):", mse_output2)

# Fit the model for Output 2
rf_regressor.fit(X_train, y_train[:, 2])
y_pred = rf_regressor.predict(X_test)
mse_output3 = mean_squared_error(y_test[:, 2], y_pred)
print("Mean Squared Error (Output 2):", mse_output3)

# %%
# Create a multi-output RandomForestRegressor
# 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[:, 2], y_pred[:, 2])
print("Mean Squared Error (Output 0):", mse_output1)
print("Mean Squared Error (Output 1):", mse_output2)
print("Mean Squared Error (Output 2):", mse_output3)

# %%
# Create a multi-output RandomForestRegressor with a scaling on the output
tt_rfr = TransformedTargetRegressor(
    regressor=RandomForestRegressor(n_estimators=n_estimators, random_state=PRNG_SEED),
    transformer=MinMaxScaler(),
)

# Fit the model
tt_rfr.fit(X_train, y_train)

# Make predictions
y_pred = tt_rfr.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[:, 2], y_pred[:, 2])
print("Mean Squared Error (Output 0):", mse_output1)
print("Mean Squared Error (Output 1):", mse_output2)
print("Mean Squared Error (Output 2):", mse_output3)
```

Labels: New Feature, help wanted