What optimization algorithms are best at traversing complicated geometries, and what trade-offs exist between different algorithms? Hi all :] so short version is just the title, and specifically as pertains to those algorithms included in R's optimx::optimx() under method = , i.e.: 'Nelder-Mead', 'BFGS', 'CG', 'L-BFGS-B', 'nlm', 'nlminb', 'spg', 'ucminf', 'newuoa', 'bobyqa', 'nmkb', 'hjkb', 'Rcgmin', and 'Rvmmin'. I have a general sense of how a few of these work, but in real-world applications some take longer to run than others and I'm not really sure when that extra time is justified or not. Presumably none of them are strictly better than any other, or else why bother including them in the package?
I'm in particular focused on which are best at traversing especially gnarly geometries full of local optima in multiple dimensions. In addition, what are good ways to improve performance & stability of optimization in this context? For example, one common trick in Bayesian phylogenetics is to raise the likelihood to some power between (0,1) to provide a more flexible proposal distribution during MC3 (i.e. 'metropolis-coupled markov chain monte carlo'). These 'heated' chains can traverse the spaces between local maxima more easily... but I'm not sure if or when that would apply using different optimization algorithms. Any ideas? I can also just re-run the optimization algorithm lotsa times and take the minimum of the runs, but that is costly.

OK, longer version of the question w/ the current motivation / background. I frequently find myself in the position of having to visualize data with some discrete property that I want to encode into the data viz using color. There are lots of discrete color palettes out there, e.g. in the packages RColorBrewer or wesanderson. Sometimes, though, I want to use a specific palette, but it has too few colors -- e.g., there are 9 colors included in "Set1" in RColorBrewer, and I need 12. So I tried to write a little script for finding new colors to add to an existing discrete color space that maximize separation from existing colors. But it looks like it tends to get stuck in local minima, judging by the fact that separate runs of optimx() return different values, each cheerfully reporting convergence. For example, here are three separate runs of my script with random initialization of color parameters from iid normals:

It’s not a ridiculously high-D optimization problem — in this application, each of my six “new colors” lives in a rescaled & logit-transformed hue-saturation-luminosity colorspace, my objective function is the harmonic mean euclidean distance to all other colors (both existing and proposed) in cartesian space (converted from cylindrical), I have a weak double exponential penalty centered on 0 to help regularize in case a color truly does occupy the boundary values so optimization doesn't chase infinity without bound, and I’m taking the minimum of multiple runs of optimx, exploring the performance of different methods included therein.
Some other notes:
In the examples above the circles are existing colors and the squares are proposed colors; the size of those points is proportional to their luminosity; and their location is given by some angle [hue] and radius [saturation] on a luminosity = 50% circle.
There's a general problem of 'non-identifiability' between proposed new colors. My attempt to solve this was found in taking my vector of colors and first exponentiating and cumulative summing all but the first of their luminosity values before inverse-logit-ing to the HSL colorspace, to constrain the first color to be the darkest and each subsequent color to be lighter than the last.
If adding just one color in 3D, a grid search could work fine, and I suppose there might be theoretical guarantees for how maximally poor performance may be for some greedy algorithm, but I don’t think finding an optimum in 18 (or however many) dimensions is too much to ask for.
I tried out other types of means (geometric, arithmetic) and other types of functions (e.g. maximizing the minimum distance to a point) but felt the harmonic to be most intuitively satisfying.
I also have few other toggles that e.g. shrink and expand the separation between points in the luminosity and saturation dimensions, since my intuition is that hue is more important for separating colors to the naked eye, and another that makes optimization prefer points that are more similar in luminosity and saturation to the existing points, since it wouldn't do to have a nice, subdued, earthy palette and get a suggestion of neon pink, or a bunch of pastels that tell you hey, dark blue would look really contrast-y here.
If I get this simpler version working properly, I might also try and weigh the visible light spectrum over the human eye's photopic luminosity function so that the angular distance might integrate over / transform according to those curves, or try to get some color-blindness weighing in there too).
Thanks for any help!

The convergence criteria just test whether the algorithm has found a local optima. As you've found, this task is not convex, so there are multiple local optima. The optimizers you're using are all local optimizers, so they will only find a globally optimal solution by chance.
When the objective function is cheap and precision isn't paramount, you can just do what you're already doing: run a local optimizer many times, and pick the best one. This is sometimes called "multi-start."
Alternatively, if the cost function is expensive, you could try one of the methods in Optimization when Cost Function Slow to Evaluate . To decide whether or not a local optimum is also a global optimum, you'll have to show that the global optimum is at least as extreme as every other local optimum. In general, this is hard. The strategies in the linked Question each have their own heuristics which are used to inform the search.
