Also, as discussed in 1:1, I think we better skip the n_jobs=1 for now and work on it as part of a second PR.

I not sure how I feel about merging a PR to main with a performance regression in the single threaded case.

It's true that on the positive side, most workloads are run with at least 2 (or 4 threads) per nodes nowadays. So maybe nobody would notice the performance regression in practice. But still...

To answer this remark in #25561 (review):

I not sure how I feel about merging a PR to main with a performance regression in the single threaded case.

As proposed in #25561 (review), we can simultaneously:

• use a feature branch to delay the inspection and the resolution
• fall-back on the previous implementation when n_jobs=1

Running this asv benchmark with this script on 1 and 16 cores confirm our intuition that we have performances increase for n_core > 1 for 32 and 64bits, on dense-dense matrices.

As discussed IRL with @jjerphan, we have performance regression for the cases sparse-sparse, dense-sparse, and sparse-dense, which could pave the way for the next PRs on this feature branch.

       before           after         ratio
     [7917117e]       [ddc464ba]
     <main>           <pull/25561/head>
-         822+/-0ms          310+/-0ms     0.38  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float32'>, 'dense', 'dense', 16)
-         823+/-0ms          265+/-0ms     0.32  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float64'>, 'dense', 'dense', 16)
       before           after         ratio
     [7917117e]       [ddc464ba]
     <main>           <pull/25561/head>
+         1.07+/-0s          3.39+/-0s     3.17  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float64'>, 'sparse', 'sparse', 1)
+         1.10+/-0s          3.38+/-0s     3.08  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float32'>, 'sparse', 'sparse', 1)
+         1.34+/-0s          3.66+/-0s     2.73  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float64'>, 'sparse', 'dense', 1)
+         185+/-0ms          505+/-0ms     2.73  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float64'>, 'sparse', 'sparse', 16)
+         193+/-0ms          506+/-0ms     2.61  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float32'>, 'sparse', 'sparse', 16)
+         1.61+/-0s          3.61+/-0s     2.24  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float32'>, 'sparse', 'dense', 1)
+         823+/-0ms          1.76+/-0s     2.14  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float64'>, 'dense', 'dense', 1)
+         257+/-0ms          540+/-0ms     2.11  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float64'>, 'sparse', 'dense', 16)
+         822+/-0ms          1.68+/-0s     2.04  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float32'>, 'dense', 'dense', 1)
+         290+/-0ms          566+/-0ms     1.96  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float32'>, 'sparse', 'dense', 16)
+         2.38+/-0s          4.52+/-0s     1.90  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float64'>, 'dense', 'sparse', 1)
+         339+/-0ms          611+/-0ms     1.80  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float64'>, 'dense', 'sparse', 16)
+         2.52+/-0s          4.53+/-0s     1.80  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float32'>, 'dense', 'sparse', 1)
+         365+/-0ms          644+/-0ms     1.76  time_pairwise_distances(10000, 10000, 10, 'manhattan',  <class 'numpy.float32'>, 'dense', 'sparse', 16)

SOME BENCHMARKS HAVE CHANGED SIGNIFICANTLY.
PERFORMANCE DECREASED.

Edit: after IRL discussion @jjerphan, we think this PR brings value because:

• Currently, on the main branch, manhattan is the only distance with sparse-sparse support. It relies on the most efficient Cython implementation, with no dispatch nor abstraction cost, which makes reconsidering the observed performance regression.
• Moreover, the remaining distance metrics don't provide sparse-sparse, sparse-dense, and dense-sparse support, nor do any other libraries in the scientific python ecosystem.
• Therefore, this PR not only improves performance in the dense-dense case but also unifies the sparse-sparse, sparse-dense, and dense-sparse support for all distance metrics.

Therefore, this PR not only improves performance in the dense-dense case but also unifies the sparse-sparse, sparse-dense, and dense-sparse support for all distance metrics.

This does not seem to always be the case in the benchmark results: there are dense-dense case with a 2x slowdown. How to you explain this?

Edit: those are the sequential cases.

Could you please re-run with the benchmarks with 2 threads?

Edit: actually, this is already visible in the second plot at the beginning of the PR.

+1 for updating the changelog as quickly as possible in PRs. It helps reviewers better understand the user-facing scope of a PR by making it explicit.

modify PairwiseDistances.is_usable_for and PairwiseDistances.valid_metrics so that the following cases aren't yet treated:
n_jobs=1

I suppose you mean the effective_n_threads measured via openmp, not the n_jobs passed by the user to a public facing API that would be meant to control the use of joblib?

Is so, +1. This means that at least in the short term we keep on using the scipy metrics implementation in the dense-dense, single threaded case and would use the new Cython infrastructure for all the other cases.

Our new benchmark is consistent with what we expected:

Note that we don't use n_jobs in this benchmark anymore, only limiting threads through threadpoolctl, which alleviates some of the overhead on main that led to the slowdown for n_jobs = 2 observed in the previous benchmark.

Hey @Micky774 ! I'm not working on it currently but I can resume this PR in the next 2 weeks if that sounds good to you?

Of course -- no rush either. I just wanted to make sure the project wasn't entirely set aside :)

Hi @Vincent-Maladiere,

#26471 (which recently has been merged in main) has introduced a separation between the interface of DistanceMetric and its implementations. This created conflicts with this PR.

Feel free to let us know if you need help resolving them. :)

On it, let's try to have something working by the end of the week :)

It seems that the merge history here got messed up @Vincent-Maladiere

@Vincent-Maladiere: if your reflog also getting in the previous state, you can find the commit of the previous state using the reflog and reset the branch to it:

# Get the latest changes from the base branch (which is up-to-date)
git checkout feature/PairwiseDistances
git pull upstream feature/PairwiseDistances

# Checkout to the previous commit
git checkout feat/pairwise_distances-pdr-backend
git reflog # find the commit
git checkout <commit>

# Force the branch of this PR back to this commit
git branch -f feat/pairwise_distances-pdr-backend

# Checkout back to this branch
git checkout feat/pairwise_distances-pdr-backend

# Update the branch, you might have conflicts to solve
git merge feature/PairwiseDistances

# After fixing conflicts, run isort on the files that you
# have modified via pre-commit.
pre-commit install
pre-commit run

# Finish the merge
git merge --continue

# Inspect the diff with the base branch
git diff feature/PairwiseDistances...

# If the diff is clear, force-push the branch 
git push -f fork feat/pairwise_distances-pdr-backend

Hey! Yes sorry about that, after rebasing on main it seems that feature/PairwiseDistance hasn't been synced for a while. @glemaitre has rebased feature/PairwiseDistance on main. So, after fixing the conflict hopefully, the large diffs will vanish.

Thanks @jjerphan for thegit branch -f trick, I'll try to fix the conflicts and simply push on the target feature branch one more time.

✔️ Linting Passed

All linting checks passed. Your pull request is in excellent shape! ☀️

Generated for commit: 54600aa

I'm running the benchmark again on {Dense, CSR} x {Dense, CSR} -- {manhattan, euclidean} distance -- {1, 2, 4, 6} threads, to make sure our conclusion is the same.

Blue is main, orange is this branch.

We obtain a slight speed-up for Euclidean distances in the Dense x CSR case (why?) and no difference otherwise, which matches our expectations, since we're not using PairwiseDistances for the Euclidean distance yet.

However, the Manhattan distance is significantly worse when using Dense x CSR or CSR x CSR, and the n_threads = 1 case with Dense Dense.

The CSR implementation we use are:

main: metrics._pairwise_fast._sparse_manhattan

branch: metrics._dist_metrics.ManhattanDistance.dist_csr

The main difference between these implementations is the parallelization happening within the _sparse_manhattan function, whereas parallelization happens outside the dist_csr in this branch, in _parallel_on_X.

Therefore, this branch makes more function call to the distance routine and this might explain the decrease.

Do you have any other ideas to explain the performance decrease? WDYT?

By the way, our tests are failing because of some connection error to conda forge.

@Vincent-Maladiere could you generate a speedscope-compatible profile report for both methods? It'll help us spot any less-than-obvious costs, as well as help us confirm valid assumptions. You can see this comment for instructions: #26316 (comment)

Okay, let's profile this. TBH, I was hoping that you or @jjerphan knew what might have impaired the benchmark so that I don't need to run Pyspy and speedscope 😅

Okay after some initial profiling, it looks like we're running into the same problems we had with the refactor attempts for the entire backend. We forfeit a lot of performance due to indirection and calls with this approach. I didn't expect it to be this dramatic, but alas.

When SparseSparseDataset64.dist is called, despite it being a rather thin wrapper around self.distance_metric.dist_csr, around 50% of the time is spent purely pushing the arguments to the stack, while the rest is the actual computation (a significant portion of which is also taken up by loading the arguments). Overall, only ~32% of time is actually spent on core computation.

This can be partially mitigated by only passing memoryviews by reference rather than value (dropping peripheral information like shape) wherever call overhead is a significant factor (e.g. in tight loops). Specifically here, it would mean changing DistanceMetric.{r}dist_csr to accept a pointer to self.{X, Y}_indices rather than the memoryview itself (#26765). This helps bridge the gap. The remaining performance loss is due to the inability to inline the computations.

I'm not sure that there is a way around this without copy-and-pasting entire swathes of the implementations. Options like code-generation and templating could solve this but at the cost of build time and binary size, along with the obvious maintenance burden.

Abstraction is directly in contest with performance here.

With that being said, this is really only a problem for manhattan distance because we have a cython specialization of it right now. We could keep using the custom specialization a special case (solely to avoid performance regression) and defer to the PairwiseDistance backend for all other cases.

@Vincent-Maladiere would you like to to keep working on this PR, or would you prefer I take over efforts here and let you focus on other things :)

Superseded by #26983

Please feel free to leave comments/suggestions on the new PR. Thank you @Vincent-Maladiere for all the work you've done so far!