[NAThompson]

Deterministic box-constrained global optimizer with the paper's local search; a small QP solver and line search replace MINQ and gls. Closes #1577

Branin function after the last sweep: boxes, base points colored by f, and the three global minima found:

Relative error vs function calls on the paper's Dixon–Szegő test set (Table II), default parameters:

Sweep-by-sweep playback on Branin and the six-hump camel, for H5Web: mcs.h5.zip

Against the library's other global optimizers: 18 problems, 25 seeds for the stochastic ones, a budget of 2000 calls per dimension, solved at f − f* ≤ max(1, |f*|)/10⁶:

All 20 problems, including rotated Rastrigin and Ackley where MCS fails: bench.h5.zip

🤖 Generated with Claude Code

[NAThompson]

@dschmitz89 : Please test at your convenience-should be quick to get it working in python for you in the nanobind layer-if it's not . . . that's another bug!

I have recently been thinking about how to validate code now that nobody reads it. I have been working on an H5Web extension which will eventually get upstreamed-but for now see here

If you run:

code --install-extension h5web-fork-deadbeef.vsix --force
code mcs.h5

there should be many compelling examples of correctness there.

[dschmitz89]

Wow, thanks again @NAThompson , you and Claude make a good team.

I will try to play around with this but will need a few days. Two questions:

• How does MCS compare to your other optimizers? Do you have benchmarks? I would expect that in terms of computation time the others will likely be be better as they do not perform linear optimizations but hopefully MCS will be able to solve more problems given a fixed function evaluation budget.
• For C++ users, it might be natural to have Eigen. Can the linear algebra here use Eigen if available and fall back to the hand woven implementation if not?

[NAThompson]

Benchmarks:

A typical run, Hartmann 6 (median of 25 seeds for the stochastic solvers):

I used 18 problems (the Dixon–Szegő set plus Rastrigin, Ackley and Rosenbrock in 2, 5 and 10 dimensions), 25 seeds for the stochastic solvers, and a budget of 2000·n calls. MCS solves 89% of the runs, against 74% for jSO, 59% for DE and 48% for CMA-ES, usually in 10–60× fewer calls. Its own overhead is about 0.5 µs per call, similar to DE and CMA-ES, so it's also fastest in total time when it solves the problem. Where it falls down: Rastrigin and Ackley in 10 dimensions, and rotated (non-separable) Rastrigin and Ackley in 5 dimensions, which jSO and CMA-ES handle better. That's expected for a method that splits along coordinates. Per-problem convergence for all 20 is in bench.h5.zip.

Eigen: the linear algebra is ≤ 0.6% of the time on most problems, and 17% of 0.11 µs per call in the worst case (Rosenbrock 10), so Eigen wouldn't help. Most of the overhead is allocation, which I can tighten later. Keeping it Eigen-free is also what lets SciPy use it; CMA-ES needs Eigen and so can't be used there.

[NAThompson]

@dschmitz89 : TLDR: This method is incredible. Thanks for the pointer!

N.B.: I am probably going to merge since this turned up some bugs in my older optimization codes (like CMA-ES) and I don't want merge conflicts. But feel free to add commentary at any time and I'll try to address it.

[dschmitz89]

Those benchmark results are awesome. One last question: the paper mentioned using tree data structures to keep RAM usage constant. Was that done here? I am asking because for example the legacy implementation of DIRECT we use in SciPy does not and indeed requires an insane amount of memory.

[NAThompson]

@dschmitz89 Good catch: no. The merged version stores each box's base point, opposite corner, split counts and model data explicitly (seven n-vectors per box) and keeps all boxes, so memory is O(calls · n). With the default budget of 1000·n calls that's 28 MB at n = 10, 0.5 GB at n = 50 and 1.1 GB at n = 100. Section 4 of the paper avoids this by storing only each box's history (parent, split index and value, child label) and recovering the rest by walking the parent chain, which makes storage proportional to the number of function evaluations. It's not constant, but it is O(calls) rather than O(calls · n). I'll implement that as a follow-up PR, which should also cut most of the solver overhead, since allocation dominated it in the profile.