Enable parallel sklearn.feature_selection.mutual_info_regression

### Describe the workflow you want to enable

I can raise the PR if someone is willing to review and potentially merge.
```
from sklearn.feature_selection import mutual_info_regression
mutual_info = mutual_info_regression(X, y, n_jobs = -1)
```

### Describe your proposed solution

In: https://github.com/scikit-learn/scikit-learn/blob/0ab36990c0a7d02052a3de7b726aa425ee14950f/sklearn/feature_selection/_mutual_info.py#L304
Change
```
    mi = [
        _compute_mi(x, y, discrete_feature, discrete_target, n_neighbors)
        for x, discrete_feature in zip(_iterate_columns(X), discrete_mask)
    ]
```
to
```
from joblib import Parallel, delayed
def process_column(x, discrete_feature):
    return _compute_mi(x, y, discrete_feature, discrete_target, n_neighbors)

mi = Parallel(n_jobs=n_jobs)(delayed(process_column)(x, discrete_feature) 
                                       for x, discrete_feature in zip(_iterate_columns(X), discrete_mask))
```

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

None

### Additional context

Enable the user to choose multicore or single core.

COMMENT:
@glemaitre I found this issue while looking for ways to accelerate `mutual_info_regression()`. It's very slow currently, and only one CPU core is used - all the others are sitting idle while my code loops very slowly from 1 to `n_features`.

Here's my benchmark. Perhaps not the benchmark you're looking for (the parallelism here is external to the function itself), but a benchmark showing how much faster this function can be when it's running in a simple parallel context. The dataset size used here is similar to an actual dataset I was looking at, which prompted my search.

```python
from sklearn.datasets import make_sparse_uncorrelated
from sklearn.feature_selection import SelectKBest, r_regression, mutual_info_regression
from sklearn.linear_model import LinearRegression
from joblib import Parallel, delayed


def objective(X, y, score_func, num_features):
    skb = SelectKBest(score_func=score_func, k=num_features)
    X_small = skb.fit_transform(X, y)
    # save X_small somewhere
    # ...
    lr = LinearRegression(n_jobs=1)
    lr.fit(X_small, y)
    return num_features, lr.score(X_small, y)


n_features = 100
n_samples = int(1e4)
X, y = make_sparse_uncorrelated(random_state=0, n_features=n_features, n_samples=n_samples)
```

Run `mutual_info_regression()` in parallel:

```python
%%time
_ = Parallel(n_jobs=-1)(
    delayed(objective)(X, y, mutual_info_regression, num_features) for num_features in range(1, n_features + 1)
)
best_res_mi_multi = dict(_)
```

```
CPU times: user 74.2 ms, sys: 111 ms, total: 186 ms
Wall time: 38 s
```

All CPU cores are used all the time.

---

Run `mutual_info_regression()` single-process:

```python
%%time
_ = Parallel(n_jobs=1)(
    delayed(objective)(X, y, mutual_info_regression, num_features) for num_features in range(1, n_features + 1)
)
best_res_mi_single = dict(_)
```

```
CPU times: user 5min 53s, sys: 1min 39s, total: 7min 33s
Wall time: 5min 31s
```

It's mostly just one CPU core that is used at any given time.

---

Comparison with `r_regression()` single-process (which may not be fair, but it's provided for context):

```python
%%time
_ = Parallel(n_jobs=1)(
    delayed(objective)(X, y, r_regression, num_features) for num_features in range(1, n_features + 1)
)
best_res_r2_single = dict(_)
```

```
CPU times: user 6.9 s, sys: 21.6 s, total: 28.5 s
Wall time: 1.79 s
```

---

My CPU:

```
$ cat /proc/cpuinfo | grep "model name" | head -n 1
model name      : AMD Ryzen 7 5800X3D 8-Core Processor
```

Labels: Needs Decision, Needs Benchmarks