Rust implementation of the SHGO algorithm (with parallel capabilities)

[jpswensen]

I have to preface this that 99% of this conversion was done with agentic AI coding. At a couple of junctures, I had tried to go in an do this by hand (with little AI help) as I was learning Rust, but finally caved and just went the route of guiding Claude Opus 4.6 to perform the conversion.

https://github.com/jpswensen/shgo-rs-native

Here was my process:

1. I guided it through the process of incrementally building up the algorithm from vertex and complex and had it write tests all along the way to make sure it gave 1-to-1 the same results as the Python implementation.
2. I guided it through the process of implementing the algorithm itself, again making sure there was 1-to-1 exact same results between the two implementations. For the delaunay triangulation, I had it use qhull. For the local optimization, I have it using the NLOPT Rust crate.
3. I then had it spend A LOT of time working on creating a modification where it parallelized the running of the various sets of points that were added to the candidate points. Parallelization is one of the areas where Rust really shines. As part of this I gave it a couple of explicit rules:
   a. If we turned off parallelization, it had to act exactly like the original Python SHGO implementation
   b. When parallelization is on, it should finish faster, but arrive at the exact same solution.
   I am pretty sure all of these were met. Sometimes, I think it does a few extra evaluations of some candidate points that might have been culled earlier in the non-parallel version, but in cases where you have 4-16 cores, these few extra evaluations of the cost function don't destroy the performance benefits of parallelization. I built in a whole suite of benchmarks for comparing the parallel and non-parallel and gave some rules of thumb about when it helps.
4. I had it build a pretty extensive test suite with a bunch of different standard functions (both hard and easy), a bunch of different dimensions, etc.
5. The initial pass had bounds capabilities, but I wanted to add linear and non-linear constraints capabilities also. I must admit that this is the point at which I have less confidence about whether there are cases where the core of SHGO's formal guarantees is fundamentally broken with the constraints. I don't think it is, because it works across many tests and gives the same results as the Python algorithm, but there may be corner cases that this completely breaks on. I think that from the original paper/code that the linear/non-linear constraints only applied to the local optimizer. This change eliminates sampled vertices that violate the bounds/constraints and also passed the constraints along to the local optimizer. I think that the only NLOPT algorithm that supports linear/non-linear constraints is COBYLA, so it is forced when constraints are added.
6. The last thing I wanted to implement was an alternative to the Delaunay triangulation. For my particular problem, once my parameter count got up to 8-10 and I wanted the number of sobol sampling points above about 64 (I think they have to be powers of 2), then it got really slow and took a lot of memory to do the Delaunay triangulation and keep track of the entire set of used and rejected candidate points. So, I dug a little bit about how some other simplex/complex optimization algorithms approach this and found that many use some nearest neighbor algorithms. So, I ended up integrating existing Rust libraries that implement a k-NN, Hierarchical Navigable Small World (HNSW, an efficient approximate nearest neighbor algorithm), and Scalable Nearest Neighbors (another approximate nearest neighbor algorithm from Google Research). I have a table below of some of the performance comparisons for the different approaches. Again, I have no idea whether this completely nukes the theoretical guarantees of SHGO, but what I can say is that on all the tests in the table below, they all arrived at the same solution as the original SHGO implementation.
7. I also got a Rust FFI wrapper so that it generates a C/C++ library. I suppose I should go make an FFI wrapper for Python also (which is made dead simple by Rust with their PyO3 crate).
8. All along the way, I ensured that if all these new features were turned off, that it gave a bit-for-bit faithful solution to the original Python implementation.

I don't even know if this is the right place to post this as an "issue", as here is the repository for the original Python implementation by the SHGO paper's authors. But, I thought it might be useful to show a couple of augmentations which are working great for me, but that I don't have a high level of confidence that it will work in ALL scenarios. IF, and that is a big if, my changes are not fundamentally algorithm breaking, I think it adds a couple of really useful capabilities in the parallelism and the Delaunay alternatives. I can tell you that the parallelism really does slam my 22 core PC at almost 80%+ utilization on all processors and does finish much faster than when running without parallelism (again see the parallel performance benchmarks). For one problem, I was able to drive my 10D problem up to n=512 and iters=8 when using the k-NN sampling method. With Delaunay, n=512 and iters=1 would have taken 24 hours just do the first sampling steps.

Any feedback would be appreciated.

Appendix:

Table of the comparative performance results from the Delaunay triangulation method and proposed alternatives.

Case	Delaunay	k-NN	HNSW	ScaNN
Rastrigin 2D	1.6 ms	0.7 ms	2.2 ms	4.4 ms
Rastrigin 3D	2.6 ms	0.8 ms	2.4 ms	4.3 ms
Ackley 5D	91.5 ms	1.7 ms	4.3 ms	16.5 ms
Rastrigin 5D	262 ms	5.6 ms	9.3 ms	65 ms
Ackley 7D	4284 ms	2.5 ms	5.6 ms	18 ms
Rastrigin 7D	17167 ms	7.5 ms	11 ms	61 ms

Table of performance comparisons for parallel vs serial

Test Case	Serial (s)	Parallel (s)	Speedup
rastrigin_5d_sobol_256	155.5	134.5	1.16x
rosenbrock_5d_sobol_256	207.8	138.4	1.50x
ackley_5d_sobol_256	135.9	123.1	1.10x
levy_5d_sobol_256	180.5	131.2	1.37x
levy_5d_sobol_512	477.6	381.3	1.25x
eggholder_2d_sobol_512	7.4	5.6	1.31x
rastrigin_2d_expensive	76.7	8.6	8.91x
rastrigin_5d_expensive	767.7	191.8	4.00x
ackley_3d_very_expensive	544.2	98.8	5.51x
eggholder_2d_sobol_512 (run 2)	7.3	5.5	1.34x

[Stefan-Endres]

Hi @jpswensen,

Thank you this is really great work! I look forward to playing around with your implementation myself.

Some questions since we have this issue open, it's maybe ideal to set up a workflow already:

1. Do you intend to maintain this repository?
2. If yes to 1, would you prefer doing future updates/features from the python version yourself, or prefer that we open pull requests to your repository? I would be happy to do the latter, ideally with help from @alchemyst since I don't have a lot of experience in Rust.

---

Other comments:

I must admit that this is the point at which I have less confidence about whether there are cases where the core of SHGO's formal guarantees is fundamentally broken with the constraints. I don't think it is, because it works across many tests and gives the same results as the Python algorithm, but there may be corner cases that this completely breaks on.

It is definitely the case that SHGO will not work on many test problems as expected. The reason is because the theory from the paper is not the same as the actual computational implementation: in particular we rely on tracking the chain algebras in the theory for things like detecting non-convex valleys (Rosenbrock function-like domains), but we not been able to implement an efficient computation of this in practice (we do still plan on implementing this).

I think that from the original paper/code that the linear/non-linear constraints only applied to the local optimizer.

Yes, intended behaviour:

1. Global simplicial complex detects Sperner minimizer domains.
2. --> constructs linear constraints (which restricts the local search domain to the detected Sperner minimizer) to pass to local minimizer.
3. --> both linear (new local and global) and non-linear constraints are passed to the local solver.

In practice (current default behaviour in implementation): linear constraints are constructed as box contraints, not the actual rigorous Sperner domain, this was done to save computation costs in higher dimensions.

The non-linear constraints are also used in the global exploration step, but currently only to eliminate sampling points that violate the constraints.

The last thing I wanted to implement was an alternative to the Delaunay triangulation. For my particular problem, once my parameter count got up to 8-10 and I wanted the number of sobol sampling points above about 64 (I think they have to be powers of 2), then it got really slow and took a lot of memory to do the Delaunay triangulation and keep track of the entire set of used and rejected candidate points. So, I dug a little bit about how some other simplex/complex optimization algorithms approach this and found that many use some nearest neighbor algorithms. So, I ended up integrating existing Rust libraries that implement a k-NN, Hierarchical Navigable Small World (HNSW, an efficient approximate nearest neighbor algorithm), and Scalable Nearest Neighbors (another approximate nearest neighbor algorithm from Google Research). I have a table below of some of the performance comparisons for the different approaches. Again, I have no idea whether this completely nukes the theoretical guarantees of SHGO, but what I can say is that on all the tests in the table below, they all arrived at the same solution as the original SHGO implementation.

Yes, this was our experience with Delaunay too, hence the development of hyperct for the default global sampling and triangulation (ensures the guarantees with the above caveats (no chain algebra tracking on valley functions)). This is still an unsolved problem in general so I think it's fine to keep exploring options here.

---

Please let me know if you have any questions and/or there's anything we can help you with.