diff --git a/src/Subsampling/benchmark/choose_n_farthest_points.cpp b/src/Subsampling/benchmark/choose_n_farthest_points.cpp
index c9c453284c..262f8d6f8b 100644
--- a/src/Subsampling/benchmark/choose_n_farthest_points.cpp
+++ b/src/Subsampling/benchmark/choose_n_farthest_points.cpp
@@ -76,7 +76,16 @@ int main(int argc, char**argv) {
             << " vs metric "              << std::chrono::duration_cast<std::chrono::milliseconds>((time_stop2 - time_start2)).count()
             << "  (Boost version " << BOOST_VERSION << ")\n";
   if (dists != dists2 || results != results2) {
+    // Note that with many points, it often happens that 2 points with the same distance are swapped in the output.
     std::cerr << "Results differ\n";
+#ifdef LOG_DIFF
+    std::ofstream log_gen("log_gen");
+    std::ofstream log_met("log_met");
+    for(std::size_t i = 0; i < results.size(); ++i){
+      log_gen << dists2[i] << '\t' << results2[i] << '\n';
+      log_met << dists [i] << '\t' << results [i] << '\n';
+    }
+#endif
     return -1;
   }
 #endif
diff --git a/src/Subsampling/include/gudhi/choose_n_farthest_points.h b/src/Subsampling/include/gudhi/choose_n_farthest_points.h
index 8b57a23c70..eaf91e676d 100644
--- a/src/Subsampling/include/gudhi/choose_n_farthest_points.h
+++ b/src/Subsampling/include/gudhi/choose_n_farthest_points.h
@@ -320,7 +320,7 @@ void choose_n_farthest_points_metric(Distance dist_,
     auto handle_neighbor_neighbors = [&](std::size_t ngb)
     {
         auto& ngb_info = landmarks[ngb];
-        std::remove_if(ngb_info.neighbors.begin(), ngb_info.neighbors.end(), [&](auto near_){
+        auto it = std::remove_if(ngb_info.neighbors.begin(), ngb_info.neighbors.end(), [&](auto near_){
             std::size_t near = near_.first;
             FT d = near_.second;
             // Conservative 3 * radius: we could use the old radii of ngb and near, but not the new ones.
@@ -328,6 +328,7 @@ void choose_n_farthest_points_metric(Distance dist_,
             // Here it is safe to use the new radii.
             return d >= max_dist(ngb_info.radius, landmarks[near].radius);
             });
+        ngb_info.neighbors.erase(it, ngb_info.neighbors.end());
     };
     // First update the Voronoi diagram, so we can compute all the updated
     // radii before pruning neighbor lists. The main drawback is that we have