ENH: use the sparse-sparse backend for computing pairwise distance

First reported in: https://github.com/scikit-learn-contrib/imbalanced-learn/issues/1056

We have a regression in `kneighbors` with sparse matrix from 1.1.X to 1.3.X.
A code sample to reproduce:

```python
# %%
import sklearn
sklearn.__version__

# %%
import numpy as np
from scipy import sparse
from sklearn.neighbors import KNeighborsRegressor

n_samples, n_features = 1_000, 10_000
X = sparse.random(n_samples, n_features, density=0.01, format="csr", random_state=0)
rng = np.random.default_rng(0)
y = rng.integers(0, 2, size=n_samples)
knn = KNeighborsRegressor(n_neighbors=5).fit(X, y)

# %%
%%timeit
knn.kneighbors(X, return_distance=False)
```

#### 1.1.X

```
21.5 ms ± 217 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
```

#### `main`

```
1.16 s ± 9.87 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```

#### Small benchmark

![image](https://github.com/scikit-learn/scikit-learn/assets/7454015/5abe644f-5c05-4124-ab43-dd0a1cdd3e58)

<details>

```python
# %%
import sklearn
sklearn.__version__

# %%
import time
from collections import defaultdict
import numpy as np
from scipy import sparse
from sklearn.neighbors import KNeighborsRegressor

n_samples, n_features = 1_000, [500, 1_000, 2_500, 5_000, 7_500, 10_000, 25_000, 50_000]

results = defaultdict(list)
for nf in n_features:
    X = sparse.random(n_samples, nf, density=0.01, format="csr", random_state=0)
    rng = np.random.default_rng(0)
    y = rng.integers(0, 2, size=n_samples)
    knn = KNeighborsRegressor(n_neighbors=5).fit(X, y)
    start = time.time()
    knn.kneighbors(X, return_distance=False)
    elapsed_time = time.time() - start
    results["version"].append(sklearn.__version__)
    results["n_features"].append(nf)
    results["elapsed_time"].append(elapsed_time)

# %%
import pandas as pd

results = pd.DataFrame(results)
results.to_csv(f"bench_{sklearn.__version__}.csv", index=False)
```

```python
# %%
import pandas as pd

results_main = pd.read_csv(f"bench_1.5.dev0.csv")
results_1_1 = pd.read_csv(f"bench_1.1.3.csv")
results = pd.concat([results_main, results_1_1], axis=0)
results

# %%
import seaborn as sns

sns.set_context("talk")

# %%
sns.relplot(
    data=results,
    x="n_features",
    y="elapsed_time",
    hue="version",
    kind="line",
    height=6,
    aspect=1.5,
    legend=True,
)
sns.despine(offset=10, trim=True)
```

</details>

#### Debug profiling

My first investigation look like we are spending all our time in the following function:

https://github.com/scikit-learn/scikit-learn/blob/5c7e831306e0a087c2b6af6913fa5b3c402f6d67/sklearn/metrics/_pairwise_distances_reduction/_middle_term_computer.pyx.tp#L39-L69

But I'm a bit rusty with profiling native code. I need a bit more time to investigate.

COMMENT:
I made an additional benchmark to check the influence of the density of the matrix:

![image](https://github.com/scikit-learn/scikit-learn/assets/7454015/77a29ed1-87c1-421e-8405-063e9c06744b)

This is also an important variable to make the switch.

Labels: Enhancement