GSOC 2026: Improve efficiency of `OptiSim` selection method

[marco-2023]

Description

OptiSim is a distance-based method implemented in Selector. The method itself scales linearly with the number of points to select. However, the algorithm is slowed down by the setup phase, where an optimal radius is searched to obtain the required number of selected data points. In the current implementation, the associated cost makes the use of large databases impractical for this method. The primary goal of this project is to fix this limitation.

📚 Package Description and Impact

The Selector library provides methods for selecting a diverse subset of a dataset. The goal is to build representative subsets from large datasets without introducing bias, and in some cases even reducing it (e.g., when the original dataset is unbalanced). This is critical for data-driven modeling across various fields and helps mitigate issues such as class imbalance in machine learning training sets. More information about this package can be found at https://www.biorxiv.org/content/10.1101/2025.11.21.689756v1 and https://selector.qcdevs.org.

👷 What will you do?

The main focus will be to refactor the OptiSim class and likely the radius optimization utility function. The objective is to decrease the computational prefactor by several orders of magnitude, while at the same time maintaining readability and a clear mapping between the implementation and the underlying mathematics.

🏁 Expected Outcomes

1. Decrease the computational cost of OptiSim by several orders of magnitude.
2. Improve the efficiency of the radius optimization used by OptiSim.
3. Write tests to ensure correctness and numerical stability after performance optimizations.
4. Provide profiling-based examples demonstrating the performance improvements.

Required skills	Python, OOP, DevOps
Preferred skills	Experience profiling and improving the performance of Python applications, preferably in scientific computing
Project size	Small
Difficulty	Medium

🙋 Mentors

Fanwang Meng	fanwang.meng_at_queensu.ca	@FanwangM
Marco Martínez-González	mmg870630_at_gmail_dot_com	@marco-2023
Farnaz Heidar-Zadeh	farnaz.heidarzadeh_at_queensu.ca	@FarnazH

🏋️ Stretch Goal

Improve the performance of the directed sphere exclusion (DISE) algorithm also.

[dhruvDev23]

Hello @marco-2023 and @FanwangM !

I have been studying the OptiSim algorithm and started working on performance optimization. I wanted to share the initial progress with you.

Problem

The original implementation was rebuilding a KDTree every iteration to find the candidate with maximum minimum distance to selected points in the algorithm method of OptiSim class in selector/methods/distance.py.

Optimization

Instead of rebuilding the KDTree repeatedly, I now maintain a min_dists array that tracks the minimum distance from each point to all selected points. This array is updated incrementally when new points are selected.

All tests passed and this optimization maintains O(N) linear scaling while improving runtime.

Please give your feedback.

[FanwangM]

Hi @dhruvDev23. Thanks for looking into this problem. The shared plots look promising. Can you please also try:

• Collect more data points for the running performance (running time vs dataset size), by starting 10, 20, 40, 60, 80, 100, 150, 200, ..., 8000.
• Plot two different cases into one plot
• Clarify how the input dataset was generated and the shape of the input data. It would be great to see how the different input data shapes affect running time (optional).

Thanks.

[dhruvDev23]

Hi @FanwangM , Apologies for the delay. I was out a little sick yesterday but am back today.

Here is the Plot of both cases (before optimization and after optimization)

Dataset Generation Details:

• Method: np.random.rand(n_points, n_features)
• Shape: (N, 10) where N ranges from 10 to 8000
• Feature space: 10-dimensional (fixed)
• Values: Uniformly distributed in [0, 1]
• Selection size: 10% of dataset size

Note: Testing with different input data shapes could be an interesting follow-up analysis. Please let me know if you'd like me to include this as well.

Please tell if you'd like to me to proceed with a PR so you can review the changes. and Should I proceed and start finding a way to optimize the optimize_radius() function.

[FanwangM]

Please go ahead and make a PR.
Also, having the analysis of different input data shapes would be helpful. Thanks.
@dhruvDev23

[PaulWAyers]

@FanwangM I talked to @FarnazH about this, and she had some ideas. It might be worth reaching out to her about it. As I recall Optisim and DISE are both supposed to be linear scaling asymptotically (but I am not at all sure; it's been a long time since I thought through them).

[dhruvDev23]

Hi @FanwangM . I have sent the PR ( #289 ). Could you please review it and let me know if the implementation looks good.

As you asked, I have also added the analysis of different input data shapes in the PR itself.

and as @PaulWAyers suggested , I will reach out to @FarnazH to ask her about the optimization ideas. Thanks.

[KumarShresth7]

Hi @FanwangM

I’m very interested in contributing to the optimization of the OptiSim algorithm and its radius optimization utility in the Selector library.

After reviewing:

selector/methods/distance.py

selector/methods/utils.py

I’ve identified the primary computational bottlenecks responsible for severe slowdowns on large datasets and have formulated a plan to reduce the computational prefactor by several orders of magnitude.

🔍 Identified Bottlenecks & Proposed Solutions
1️⃣ Setup Phase Bottleneck (optimize_radius)
The Problem

optimize_radius currently performs a binary search for the optimal radius r by repeatedly calling:

obj.algorithm(X, upper_size)

This is done up to n_iter times.

Running the full selection process just to evaluate a single radius guess is computationally prohibitive for large datasets.

Proposed Solution

Instead of running the full OptiSim algorithm for every radius guess:

Use spatial.KDTree(X) to estimate average neighborhood density.

Perform queries on a randomized subsample to compute a highly accurate initial guess for r.

Replace full runs during search with a lightweight greedy exclusion loop that:

Only counts excluded points

Skips expensive k-candidate selection logic

This allows:

Rapid bounding of radius values

Only one final full OptiSim execution

Orders-of-magnitude reduction in setup cost

2️⃣ Inner Loop Bottleneck (OptiSim.algorithm)
The Problem

Inside the while loop, the code rebuilds a KDTree on every iteration:

new_tree = spatial.KDTree(x[selected])

As selected size S grows:

Tree construction becomes O(S log S)

Done S times → O(S² log S) just for distance lookups

This dominates runtime for large selections.

Proposed Solution

Since the candidate sublist is restricted to size k (default = 10):

Instead of rebuilding a KDTree:

Compute distances directly using:

scipy.spatial.distance.cdist (vectorized)

OR maintain an incrementally updated minimum-distance array

This transforms:

O(S² log S) → O(k · S)

Benefits:

Eliminates unnecessary tree reconstruction

Dramatically reduces inner-loop overhead

Keeps logic mathematically equivalent

I would appreciate your thoughts on:

Whether the lightweight radius approximation should:

Be tightly integrated into utils.py

OR abstracted into a separate heuristic module

Thank you for your time and consideration.
I look forward to the opportunity to contribute to the Selector package.

Best regards,
Kumar Shresth

[PaulWAyers]

@KumarShresth7 in general I would imagine that a "utility" that isn't tightly tied to the "mission" of selector would go into utils.py. The exception is if the utility is quite large, in which case a new module might be justifiable. Overall, though, the goal is to keep the details that the average user won't need to know about as separate from the main code as possible, while ensuring those details are "somewhere one expects" so they can be found as needed. We also try to keep things (even with a module) as modular as possible, so that (inevitable!) failures are localized and easy to detect.

[KumarShresth7]

Thanks for the clarification!

That makes perfect sense. I will proceed by keeping the implementation within selector/methods/utils.py.

To align with your points on modularity and "hiding details":

I will refactor the existing optimize_radius function in utils.py.

I will implement the new "density approximation" logic as a distinct, private helper function (e.g., _estimate_radius_density) within that same module.

This keeps the main OptiSim class clean and ensures the heavy lifting happens "somewhere one expects" without cluttering the user-facing API.

I’ll ensure these components are loosely coupled so that if the approximation fails or needs adjustment, it’s easy to isolate. I’m ready to start on the implementation!

[dhruvDev23]

Hi @FanwangM , @marco-2023 and @FarnazH ,

My initial attempt at a min_dists array optimization approach was not quite right for pairwise datasets. so I've refocused on refactoring the algorithm() method of OptiSim method directly to address the core bottlenecks found in the original implementation.

I did the profiling again with the profiling code provided by Marco. The results point to the query_ball_point calls as the time consuming step. and he suggested to try to refactor the algorithm method to decrease query_ball_point use or its input.

Here are the Profiling results of OptiSim class with the original implementation

1082931 function calls (1080745 primitive calls) in 29.648 seconds

   Ordered by: cumulative time
   List reduced from 148 to 30 due to restriction <30>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000   29.648   29.648 base.py:38(select)
        1    0.000    0.000   29.648   29.648 distance.py:411(select_from_cluster)
        1    0.001    0.001   29.647   29.647 utils.py:35(optimize_radius)
       11    1.436    0.131   29.621    2.693 distance.py:339(algorithm)
     2195   21.897    0.010   21.905    0.010 _kdtree.py:480(query_ball_point)
     2204    2.853    0.001    3.499    0.002 _kdtree.py:353(__init__)
     2193    0.136    0.000    2.527    0.001 _kdtree.py:363(query)
    87535    1.447    0.000    1.447    0.000 {method 'acquire' of '_thread.lock' objects}
    17507    0.044    0.000    1.374    0.000 threading.py:873(start)
    17507    0.045    0.000    0.847    0.000 threading.py:556(wait)
    17507    0.022    0.000    0.835    0.000 threading.py:1021(join)
    17507    0.017    0.000    0.800    0.000 threading.py:1059(_wait_for_tstate_lock)
    17507    0.058    0.000    0.778    0.000 threading.py:280(wait)
    15388    0.661    0.000    0.661    0.000 {method 'reduce' of 'numpy.ufunc' objects}
     8794    0.016    0.000    0.622    0.000 fromnumeric.py:69(_wrapreduction)
    17507    0.478    0.000    0.478    0.000 {built-in method _thread.start_new_thread}
     2204    0.003    0.000    0.304    0.000 fromnumeric.py:2903(amax)
     2204    0.002    0.000    0.290    0.000 fromnumeric.py:3046(amin)
     2195    0.011    0.000    0.122    0.000 core.py:6779(array)
     2195    0.028    0.000    0.111    0.000 core.py:2833(__new__)
    17507    0.037    0.000    0.106    0.000 threading.py:795(__init__)
     4400    0.009    0.000    0.081    0.000 core.py:3893(compressed)
     4400    0.069    0.000    0.069    0.000 {method 'compress' of 'numpy.ndarray' objects}
     6592    0.002    0.000    0.055    0.000 _methods.py:67(_all)
8779/6595    0.005    0.000    0.046    0.000 {method 'view' of 'numpy.ndarray' objects}
     2195    0.027    0.000    0.042    0.000 core.py:3004(__array_finalize__)
    17507    0.008    0.000    0.038    0.000 threading.py:521(__init__)
     2193    0.004    0.000    0.030    0.000 fromnumeric.py:3068(prod)
     8798    0.027    0.000    0.027    0.000 {built-in method numpy.array}
     6579    0.006    0.000    0.027    0.000 fromnumeric.py:51(_wrapfunc)

Breakdown of the most time consuming processes:

1. query_ball_point: 21.905 seconds (74% of total execution time), 2195 calls

• The algorithm calls this repeatedly inside the while loop to eliminate candidates within radius.
• At this scale, the overhead of multi-threading (workers=-1) and lock contention far exceeds the actual search cost, triggering over 17,000 thread spawns.

2. KDTree Initialization: 2.853 seconds, 2204 calls

• The original code instantiates a brand new KDTree in every single iteration of the selection loop.

Proposed Solution:
I will refactor algorithm() method of OptiSim and create two optimized paths for two types of inputs:

• Path A — Distance Matrix Input : If the input x is already a square distance matrix, the algorithm would use direct NumPy row slicing: bv[X_dist[idx] <= self.r] = 1. This would eliminate query_ball_point and KDTree usage entirely for precomputed data.

• Path B — Raw Feature Input : Instead of rebuilding a new KDTree in every loop, we can replace the inner loop search with a direct NumPy distance calculation between the candidates and the selected points. This local calculation is extremely fast and avoids the O(N^2) memory issues of a full precomputed matrix.

These changes would significantly optimize OptiSim and make it more practical for large datasets.

Please let me know if I should start working on this implementation and do profiling and benchmarks.

I would love to hear your feedback on this approach :)

[dhruvDev23]

Hi Mentors,

I have implemented the two paths approach that I proposed earlier for algorithm method of OptiSim class. I have done the benchmarks using the benchmark code provided by marco. Below are the benchmark results done using 2000 point dataset (size=25, averaged over 10 runs):

Benchmark Results:

Input Type : Pairwise Distance Matrix

• Before optimization

Mean runtime: 0.8016 s
Std dev:      0.0082 s
Min runtime:  0.7911 s
Max runtime:  0.8196 s

• After optimization

Mean runtime: 0.0090 s
Std dev:      0.0005 s
Min runtime:  0.0086 s
Max runtime:  0.0105 s

~89x Speed Up

Input Type : Raw Feature Matrix

• Before optimization

Mean runtime: 0.0296 s                                 
Std dev:      0.0017 s                                       
Min runtime:  0.0283 s                                   
Max runtime:  0.0333 s                                  

• After optimization

Mean runtime: 0.0058 s
Std dev:      0.0004 s
Min runtime:  0.0053 s
Max runtime:  0.0070 s

~5x Speed Up

Changes:

• Path A (Pairwise Distance Matrix): Eliminated KDTree and query_ball_point entirely. Exclusion uses vectorized NumPy boolean masking bv[x_dist[idx] <= self.r] = 1, and candidate evaluation uses direct matrix slicing x_dist[np.ix_(selected, sublist)].min(axis=0).

• Path B (Raw Feature Matrix): Kept the single outer KDTree for exclusion but eliminated the inner-loop KDTree rebuild. Replaced it with a direct NumPy distance calculation np.linalg.norm.

All existing tests pass with no regressions.

Please let me know if this approach and these results are reasonable. If so, I would love to move on to next step. Based on the project description and the call stack in the profiling results, my next target will be to increase efficiency of optimize_radius function in utils.py.

[dhruvDev23]

Hi @PaulWAyers , @FanwangM and @marco-2023 ,

I am currently preparing my GSoC proposal for the "Improve efficiency of OptiSim selection method" project, and I have a couple of quick questions regarding the scope:

1. Project Size Discrepancy: I noticed that the GitHub Issues section lists the project as "Small" (90 hours), whereas the QC Devs website lists it as "Medium" (175 hours). Could you please clarify which project size I should build my proposal timeline around?
2. Stretch Goals: I am also interested in the stretch goal to "Improve the performance of the directed sphere exclusion (DISE) algorithm". If this is intended to be a Medium (175-hour) project, would you recommend integrating the DISE optimization into my main proposal timeline, or should I keep it as a stretch goal to pursue only if the OptiSim work is completed ahead of schedule?

Thank you for your help!

Best regards,
Dhruv

[ryusufsadiq-spec]

Hello! I am a 300-level Mechanical Engineering student at the University of Abuja. I am very interested in this project for GSoC 2026 as it aligns with my interest in optimization and physics-based data. I would love to contribute. Are there any specific warm-up tasks I should start with?