\chisao{}: A GPU-Native Parallel Optimizer for Multimodal Black-Box Functions via Convergence-Anticonvergence Oscillation

# A GPU-Native Parallel Optimizer for Multimodal Black-Box Functions via Convergence-Anticonvergence Oscillation
Source: [https://arxiv.org/html/2606.26164](https://arxiv.org/html/2606.26164)
Ira Wolfson Department of Electronics and Electrical Engineering Braude College of Engineering, Karmiel, Israel wolfsoni@braude\.ac\.il

###### Abstract

Finding all modes of a multimodal black\-box function is a fundamental challenge in optimization, Bayesian inference, and scientific computing\. Existing approaches—basin\-hopping, CMA\-ES, multistart gradient descent—operate sequentially and cannot exploit the massive parallelism of modern GPU hardware\. We introduceChiSao\(Convergence\-Halt\-Invert\-Stick\-And\-Oscillate\), a GPU\-native population optimizer that runs an entire sample batch simultaneously and exploits a deliberate convergence\-anticonvergence oscillation cycle to escape local traps while freezing confirmed modes\. The structural move is asymmetric: samples that reach true peaks are frozen \(“stuck”\) and preserved, while the rest keep exploring via momentum\-based anti\-convergence and stochastically smoothed gradients\. Adaptive reseeding via two complementary strategies \(Repulse Monkey and Golden Rooster\) maintains population diversity throughout\. On all 42 functions of the Simon Fraser University optimization benchmark suite across dimensionsd∈\{2,4,8,16,32,64\}d\\in\\\{2,4,8,16,32,64\\\},ChiSaoachieves100%mode recovery where all CPU baselines collapse atd≥8d\\geq 8on the hardest multimodal functions, at up to34×34\\timesspeedup over basin\-hopping on functions where all methods succeed \(Michalewiczd=64d=64\) and up to39×39\\timeson unimodal functions \(Rotated Hyper\-Ellipsoidd=64d=64, pure GPU dividend\)\. All benchmarks evaluate the objective by value alone—gradients come from finite differences—so the reported speedups are a derivative\-free worst case\. Under substantial likelihood noise \(σnoise\\sigma\_\{\\mathrm\{noise\}\}up to 1\.0\), mode detection remains 100% reliable\. The algorithm is available as a standalone open\-source Python package on PyPI\.

Keywords:multimodal optimization, GPU computing, black\-box optimization, population methods, mode finding, parallel optimization

## 1Introduction

Multimodal optimization—finding all significant maxima of a functionf:ℝd→ℝf:\\mathbb\{R\}^\{d\}\\to\\mathbb\{R\}accessible only through black\-box evaluations—arises across scientific computing, probabilistic inference, and machine learning\. The challenge is not merely finding a global optimum but*cataloguing all modes*, since downstream tasks \(Bayesian model averaging, multistart refinement, mixture fitting\) require knowledge of the full mode structure\.

Classical approaches fall into two families\.*Sequential*methods such as basin\-hopping\[[27](https://arxiv.org/html/2606.26164#bib.bib1)\]and simulated annealing\[[10](https://arxiv.org/html/2606.26164#bib.bib3)\]perturb a single solution and accept or reject moves stochastically; they explore one trajectory at a time and are inherently sequential\.*Population*methods such as CMA\-ES\[[6](https://arxiv.org/html/2606.26164#bib.bib6)\]and particle swarm optimization\[[9](https://arxiv.org/html/2606.26164#bib.bib7)\]maintain a set of candidate solutions but typically converge the entire population toward a single basin, losing mode diversity\. Multistart gradient descent runs multiple independent optimization trajectories but these are embarrassingly parallel only in the trivial sense—each trajectory is independent and discovers at most one mode\.

None of these is built for the architecture of modern GPUs: thousands of arithmetic units running the same instruction on different data simultaneously\. On a GPU withP\>104P\>10^\{4\}parallel cores, a batch ofN≤PN\\leq Pgradient evaluations completes in the same wall\-clock time as a single evaluation\. Sequential methods waste this capacity entirely; population methods use it only within a single generation step\.

ChiSao111Wing Chun*Chí Sǎo*\(“Sticky Hands”\): in the martial art, sensitivity drills train simultaneous contact maintenance and redirection, precisely the explore\-and\-freeze asymmetry of the algorithm\.is a GPU\-native population optimizer with three structural moves no prior method combines\.

The first is a freeze\-and\-explore asymmetry\. A sample that reaches a true peak, passing both a gradient\-norm and a likelihood\-quality test, is frozen from the exploration phases but continues to participate in gradient ascent for peak refinement\. The rest of the population keeps moving\. Standard population methods either advance all particles or none;ChiSaoadvances only those that still have work to do\.

The second is a deliberate anti\-convergence phase\. After each batched L\-BFGS pass, unfrozen samples take momentum\-based gradient*descent*steps\. This is neither noise injection \(simulated annealing\) nor random perturbation \(basin\-hopping\); it is directed motion, the momentum term carrying samples across valleys into new basins\.

The third is stochastic smoothing, applied as a Hands Like Clouds phase in which unfrozen samples ascend on a Gaussian\-smoothed estimate of∇f\\nabla f\. Smoothing at scaleσ\\sigmaerases sub\-σ\\sigmatexture and exposes the global basin geometry\. HLC and anti\-convergence are complementary rather than redundant: HLC steers samples toward promising regions on the smoothed landscape; anti\-convergence then disperses them on the raw landscape\.

Two reseeding strategies—Repulse Monkey when many samples remain unfrozen, Golden Rooster when most are frozen—close the cycle\. The algorithm needs only batched function and gradient evaluations, maps directly onto GPU execution, and works as a drop\-in exploration module for any optimizer that supplies an initial population\.

##### Contributions\.

The contribution is fourfold: the oscillation cycle as a strategy for multimodal optimization with a formalized freeze\-and\-explore asymmetry \(Section[3](https://arxiv.org/html/2606.26164#S3)\); convergence guarantees for log\-concave targets and a coverage analysis for the multimodal case \(Section[4](https://arxiv.org/html/2606.26164#S4)\); benchmarks against Differential Evolution, Basin\-Hopping, and CMA\-ES on all 42 SFU functions acrossd∈\{2,4,8,16,32,64\}d\\in\\\{2,4,8,16,32,64\\\}\(Section[5](https://arxiv.org/html/2606.26164#S5)\); and noise\-robustness results showing 100% mode detection up toσnoise=1\.0\\sigma\_\{\\mathrm\{noise\}\}=1\.0, the signal scale itself \(Section[5\.7](https://arxiv.org/html/2606.26164#S5.SS7)\)\.

## 2Related Work

##### Basin\-hopping\.

Wales and Doye\[[27](https://arxiv.org/html/2606.26164#bib.bib1)\]introduced basin\-hopping as alternating random perturbations with local minimization, accepting moves via a Metropolis criterion\. It is effective in low dimensions and has been widely used in computational chemistry\[[28](https://arxiv.org/html/2606.26164#bib.bib2)\]\. Its fundamental limitation is sequential execution: one perturbation, one local minimization, one accept/reject decision per step\. GPU parallelization is limited to the local minimization sub\-problem, not the search strategy\.

##### Simulated annealing\.

Kirkpatrick et al\.\[[10](https://arxiv.org/html/2606.26164#bib.bib3)\]introduced temperature\-based stochastic acceptance, which provides theoretical convergence guarantees under slow cooling\[[5](https://arxiv.org/html/2606.26164#bib.bib4)\]but requires an exponentially slow schedule to avoid premature convergence in practice\. Parallel implementations\[[21](https://arxiv.org/html/2606.26164#bib.bib5)\]distribute independent chains but do not share information, losing the benefit of population diversity\.

##### Evolutionary and population methods\.

CMA\-ES\[[6](https://arxiv.org/html/2606.26164#bib.bib6)\]adapts a full covariance matrix to the population geometry, achieving state\-of\-the\-art performance on unimodal and mildly multimodal problems\. Niching extensions\[[20](https://arxiv.org/html/2606.26164#bib.bib11)\]attempt to maintain diversity but add significant complexity\. Particle swarm optimization\[[9](https://arxiv.org/html/2606.26164#bib.bib7)\]uses social attraction toward the population best, which concentrates the population at a single basin in multimodal settings unless modified with repulsion terms\[[3](https://arxiv.org/html/2606.26164#bib.bib8)\]\. Differential evolution\[[25](https://arxiv.org/html/2606.26164#bib.bib9)\]generates candidates by combining population members but shares CMA\-ES’s tendency toward mode collapse\.

##### Multistart methods\.

Running multiple independent gradient\-based optimizers from random initial points is the most common approach in practice\[[13](https://arxiv.org/html/2606.26164#bib.bib10)\]\. It parallelizes trivially but has two weaknesses: restarts that converge to the same mode waste computation, and there is no mechanism to ensure the population maintains diversity after convergence begins\.ChiSao’s deduplication and reseeding directly address both\.

##### Bayesian optimization\.

Gaussian\-process surrogate methods\[[24](https://arxiv.org/html/2606.26164#bib.bib25)\]target the same black\-box setting from a serial angle: each acquisition step refits the surrogate against all prior evaluations and then optimises the acquisition function\. The design objective is locating the single best point at minimum evaluation budget, not cataloguing modes, and the per\-step surrogate cost grows poorly with population size\.

##### GPU\-parallel optimization\.

GPU acceleration of optimization has focused primarily on stochastic gradient descent for deep learning\[[4](https://arxiv.org/html/2606.26164#bib.bib12)\], where the objective is unimodal and the bottleneck is data throughput rather than exploration\. For black\-box multimodal optimization, GPU acceleration is largely unexplored, withSalimanset al\.\[[23](https://arxiv.org/html/2606.26164#bib.bib24)\]the major exception—evolution strategies massively parallelized for reinforcement learning, but with a single\-objective design that does not catalogue modes\.Albert \[[1](https://arxiv.org/html/2606.26164#bib.bib13)\]parallelized nested sampling on GPU \(JAXNS\) but the sequential compression structure of nested sampling limits GPU utilization to within\-iteration batching\.ChiSaoappears to be the first black\-box multimodal optimizer with full GPU utilization across the search strategy itself, not only inside individual iterations\.

##### Smoothing and continuation methods\.

Gaussian smoothing of objective functions as a heuristic for escaping local optima has been studied under the name “diffusion” or “graduated non\-convexity”\[[2](https://arxiv.org/html/2606.26164#bib.bib14),[14](https://arxiv.org/html/2606.26164#bib.bib15)\]\.Nesterov and Spokoiny \[[16](https://arxiv.org/html/2606.26164#bib.bib16)\]analyzed Gaussian smoothing for convex optimization\. The Hands Like Clouds phase ofChiSaoapplies smoothing as one phase within a structured oscillation cycle, rather than as a global preprocessing step or annealing schedule\.

##### Position relative to prior work\.

ChiSaodiffers from every method above on at least one structural axis: directional anti\-convergence rather than random perturbation, stochastic smoothing within a cycle rather than a global annealing schedule, freeze\-and\-explore asymmetry rather than uniform population movement, and full\-batch GPU execution rather than within\-step parallelism\. Basin\-hopping with restarts\[[27](https://arxiv.org/html/2606.26164#bib.bib1)\]is the closest neighbour; it differs on all four\.

## 3Algorithm

### 3\.1Problem Statement

Letf:Θ→ℝf:\\Theta\\to\\mathbb\{R\}be a black\-box function on a bounded domainΘ⊂ℝd\\Theta\\subset\\mathbb\{R\}^\{d\}, accessible through evaluationsf​\(θ\)f\(\\theta\); the gradients∇f​\(θ\)\\nabla f\(\\theta\)the algorithm uses are computed by finite differences, or supplied analytically when available\. We seek the set of all*significant modes*:

ℳ∗=\{θ∗∈Θ:∇f​\(θ∗\)=0,∇2f​\(θ∗\)≺0,f​\(θ∗\)≥fmax\+log⁡δ\}\\mathcal\{M\}^\{\*\}=\\left\\\{\\theta^\{\*\}\\in\\Theta:\\nabla f\(\\theta^\{\*\}\)=0,\\;\\nabla^\{2\}f\(\\theta^\{\*\}\)\\prec 0,\\;f\(\\theta^\{\*\}\)\\geq f\_\{\\max\}\+\\log\\delta\\right\\\}\(1\)wherefmax=maxθ⁡f​\(θ\)f\_\{\\max\}=\\max\_\{\\theta\}f\(\\theta\)andδ∈\(0,1\)\\delta\\in\(0,1\)is a quality threshold \(defaultδ=0\.1\\delta=0\.1, i\.e\., modes within one decade of the global maximum\)\. The goal is not merely to find the global maximum but to identify all members ofℳ∗\\mathcal\{M\}^\{\*\}\.

ChiSaotakes as input a population\{xi\}i=1N⊂Θ\\\{x\_\{i\}\\\}\_\{i=1\}^\{N\}\\subset\\Thetaof initial candidate points \(arbitrary; provided by the caller\), and returns the estimated setℳ^∗⊆\{xi\}\\hat\{\\mathcal\{M\}\}^\{\*\}\\subseteq\\\{x\_\{i\}\\\}\.

### 3\.2The Oscillation Cycle

ChiSaoruns fornoscn\_\{\\mathrm\{osc\}\}oscillation cycles \(default 3\)\. Each cycle consists of six phases executed in fixed order\. The ordering is deliberate and is justified in Section[3\.3](https://arxiv.org/html/2606.26164#S3.SS3)\.

#### 3\.2\.1Phase 1: Convergence

All stuck masks are released \(every sample is temporarily unfrozen\) and the entire population undergoesnconvn\_\{\\mathrm\{conv\}\}steps of batched L\-BFGS\[[12](https://arxiv.org/html/2606.26164#bib.bib17),[17](https://arxiv.org/html/2606.26164#bib.bib27)\]toward local maxima offf\. AllNNsamples optimize simultaneously in a single GPU batch:

xi←L​\-​BFGS​\(xi,nconv\)∀i=1,…,N\(GPU\-parallel\)x\_\{i\}\\leftarrow\\mathrm\{L\\text\{\-\}BFGS\}\(x\_\{i\},n\_\{\\mathrm\{conv\}\}\)\\quad\\forall i=1,\\ldots,N\\quad\\text\{\(GPU\-parallel\)\}\(2\)The L\-BFGS memory parameter ism=10m=10by default\. Iteration count adapts to dimension:nconv=max⁡\(10,3​log2⁡d\)n\_\{\\mathrm\{conv\}\}=\\max\(10,3\\log\_\{2\}d\)\.

#### 3\.2\.2Phase 2: Stick Detection

After convergence, samples that have reached true local maxima are marked as*stuck*:

stucki←\(‖∇f​\(xi\)‖∞<ϵgrad\)∧\(f​\(xi\)≥fmax\+log⁡δ\)\\mathrm\{stuck\}\_\{i\}\\leftarrow\\left\(\\\|\\nabla f\(x\_\{i\}\)\\\|\_\{\\infty\}<\\epsilon\_\{\\mathrm\{grad\}\}\\right\)\\wedge\\left\(f\(x\_\{i\}\)\\geq f\_\{\\max\}\+\\log\\delta\\right\)\(3\)The gradient thresholdϵgrad\\epsilon\_\{\\mathrm\{grad\}\}\(default10−610^\{\-6\}\) ensures the sample has reached a stationary point\. The likelihood thresholdlog⁡δ\\log\\delta\(defaultlog⁡0\.1\\log 0\.1\) prevents sticking at low\-quality local maxima that fall outsideℳ∗\\mathcal\{M\}^\{\*\}\. Stuck samples are excluded from exploration phases but continue to participate in gradient ascent for peak refinement\.

#### 3\.2\.3Phase 3: Deduplication

Stuck samples are deduplicated using theL∞L\_\{\\infty\}metric \(see Appendix[A](https://arxiv.org/html/2606.26164#A1)\):

remove​xj​if​∃xi≠xj:‖xi−xj‖∞<ϵdup,f​\(xi\)≥f​\(xj\)\\text\{remove \}x\_\{j\}\\text\{ if \}\\exists\\,x\_\{i\}\\neq x\_\{j\}:\\\|x\_\{i\}\-x\_\{j\}\\\|\_\{\\infty\}<\\epsilon\_\{\\mathrm\{dup\}\},\\;f\(x\_\{i\}\)\\geq f\(x\_\{j\}\)\(4\)The sample with the higher function value survives; its inverse Hessian estimate is preserved for width estimation\. The number of removed duplicatesKlostK\_\{\\mathrm\{lost\}\}is recorded\.

#### 3\.2\.4Phase 4: Reseeding \(Repulse Monkey / Golden Rooster\)

IfKlost\>0K\_\{\\mathrm\{lost\}\}\>0and this is not the last oscillation, the population is maintained at sizeNNby reseeding\. Two strategies are used depending on exhaustion state\.

##### Repulse Monkey\.

When≥5\\geq 5unfrozen samples remain, new samples are generated by shooting rays from unfrozen samples in random directions, sampling uniformly along those rays withinΘ\\Theta\. This disperses new candidates away from known peaks\.

##### Golden Rooster\.

When<5<5unfrozen samples remain \(near exhaustion\), new samples are generated from confirmed peaks using orthonormal ray directions obtained via QR decomposition of a random matrix\. This exploits the GPU’s parallel capacity to systematically probe orthogonal directions from known peaks\. Golden Rooster adaptively disables itself if successive generations discover no new peaks, preventing wasted computation\.

All reseeded samples are marked as unstuck to enable exploration\.

#### 3\.2\.5Phase 5: Hands Like Clouds

If this is not the last oscillation, unfrozen samples takencloudn\_\{\\mathrm\{cloud\}\}gradient*ascent*steps using a stochastically smoothed gradient:

∇fσ​\(θ\)≈1K​∑k=1K∇f​\(θ\+σ​zk\),zk∼𝒩​\(0,Id\)\\nabla f\_\{\\sigma\}\(\\theta\)\\approx\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\\nabla f\(\\theta\+\\sigma z\_\{k\}\),\\quad z\_\{k\}\\sim\\mathcal\{N\}\(0,I\_\{d\}\)\(5\)Smoothing at scaleσ\\sigmablurs features smaller thanσ\\sigmawhile preserving global basin structure\. For multi\-scale problems \(Rastrigin, Ackley\), this allows samples to ascend toward global basins through local ripples\. The smoothing scaleσ\\sigmais auto\-estimated from sample spread:σ=c⋅std​\(\{xi:¬stucki\}\)\\sigma=c\\cdot\\mathrm\{std\}\(\\\{x\_\{i\}:\\neg\\mathrm\{stuck\}\_\{i\}\\\}\)withc=0\.1c=0\.1by default, or supplied by the caller\. Stuck samples are excluded from this phase\.

#### 3\.2\.6Phase 6: Anti\-Convergence

Unfrozen samples receivenantin\_\{\\mathrm\{anti\}\}momentum\-based gradient*descent*steps:

vt\+1=μ⋅vt−α⋅∇f​\(xt\),xt\+1=xt\+vt\+1v\_\{t\+1\}=\\mu\\cdot v\_\{t\}\-\\alpha\\cdot\\nabla f\(x\_\{t\}\),\\qquad x\_\{t\+1\}=x\_\{t\}\+v\_\{t\+1\}\(6\)with momentumμ=0\.9\\mu=0\.9and step sizeα\\alpha\(auto\-tuned from scale estimates or provided by caller\)\. Descent in the gradient direction moves samples away from current peaks; the momentum term carries samples across valleys and into new basins\. Stuck samples are excluded from this phase\.

Phases 4–6 are skipped on the last oscillation \(k=nosck=n\_\{\\mathrm\{osc\}\}\), so that the final convergence phase yields clean output peaks without further perturbation\.

### 3\.3Design Rationale

The phase ordering is not arbitrary\. Deduplication \(Phase 3\) follows convergence \(Phase 1\) because samples are maximally clustered at confirmed peaks after convergence, making duplicate detection most effective\. Reseeding \(Phase 4\) precedes smoothing \(Phase 5\) so that newly injected samples benefit from the smoothed\-gradient phase\. Smoothing \(Phase 5\) precedes anti\-convergence \(Phase 6\) to nudge samples toward global basins before momentum\-based exploration begins; reversing this order would apply smoothing to already\-dispersed samples, wasting its effect\. Anti\-convergence \(Phase 6\) is the last exploration step so that samples enter the next convergence phase from maximally diverse positions\.

The release of all stuck masks at the start of Phase 1 \(“Release all”\) is also deliberate: it allows samples that were previously confirmed but may have been perturbed by reseeding to re\-converge from their current positions, rather than remaining locked at stale peak estimates\.

### 3\.4Complete Algorithm

Algorithm 1ChiSao: Convergence\-Anticonvergence Oscillation Optimizer1:Initial population

\{xi\}i=1N⊂Θ\\\{x\_\{i\}\\\}\_\{i=1\}^\{N\}\\subset\\Theta, function

ff\(gradients

∇f\\nabla fby finite differences, or analytic if supplied\)

2:Parameters:

nosc,nconv,nanti,ncloud,ϵgrad,ϵdup,δ,μ,α,σn\_\{\\mathrm\{osc\}\},n\_\{\\mathrm\{conv\}\},n\_\{\\mathrm\{anti\}\},n\_\{\\mathrm\{cloud\}\},\\epsilon\_\{\\mathrm\{grad\}\},\\epsilon\_\{\\mathrm\{dup\}\},\\delta,\\mu,\\alpha,\\sigma
3:Estimated mode set

ℳ^∗\\hat\{\\mathcal\{M\}\}^\{\*\}
4:

stuck←𝐅𝐚𝐥𝐬𝐞N\\mathrm\{stuck\}\\leftarrow\\mathbf\{False\}^\{N\}
5:for

k=1,…,nosck=1,\\ldots,n\_\{\\mathrm\{osc\}\}do

6:

stuck←𝐅𝐚𝐥𝐬𝐞N\\mathrm\{stuck\}\\leftarrow\\mathbf\{False\}^\{N\}⊳\\trianglerightPhase 1: release all, then converge

7:

x,∇f​\(x\)←L​\-​BFGS​\_​batch​\(x,nconv\)x,\\nabla f\(x\)\\leftarrow\\mathrm\{L\\text\{\-\}BFGS\\\_batch\}\(x,\\,n\_\{\\mathrm\{conv\}\}\)
8:

fmax←maxi⁡f​\(xi\)f\_\{\\max\}\\leftarrow\\max\_\{i\}f\(x\_\{i\}\)
9:

stucki←\(‖∇f​\(xi\)‖∞<ϵgrad\)∧\(f​\(xi\)≥fmax\+log⁡δ\)\\mathrm\{stuck\}\_\{i\}\\leftarrow\\bigl\(\\\|\\nabla f\(x\_\{i\}\)\\\|\_\{\\infty\}<\\epsilon\_\{\\mathrm\{grad\}\}\\bigr\)\\wedge\\bigl\(f\(x\_\{i\}\)\\geq f\_\{\\max\}\+\\log\\delta\\bigr\)⊳\\trianglerightPhase 2

10:

x,Klost←deduplicate​\(x,stuck,ϵdup\)x,\\,K\_\{\\mathrm\{lost\}\}\\leftarrow\\mathrm\{deduplicate\}\(x,\\,\\mathrm\{stuck\},\\,\\epsilon\_\{\\mathrm\{dup\}\}\)⊳\\trianglerightPhase 3

11:if

k<nosck<n\_\{\\mathrm\{osc\}\}then

12:

x←reseed​\(x,stuck,Klost,Θ\)x\\leftarrow\\mathrm\{reseed\}\(x,\\,\\mathrm\{stuck\},\\,K\_\{\\mathrm\{lost\}\},\\,\\Theta\)⊳\\trianglerightPhase 4: Repulse Monkey / Golden Rooster

13:

x​\[¬stuck\]←HLC​\(x​\[¬stuck\],ncloud,σ\)x\[\\neg\\mathrm\{stuck\}\]\\leftarrow\\mathrm\{HLC\}\(x\[\\neg\\mathrm\{stuck\}\],\\,n\_\{\\mathrm\{cloud\}\},\\,\\sigma\)⊳\\trianglerightPhase 5: Hands Like Clouds

14:

x​\[¬stuck\]←anti​\_​converge​\(x​\[¬stuck\],nanti,μ,α\)x\[\\neg\\mathrm\{stuck\}\]\\leftarrow\\mathrm\{anti\\\_converge\}\(x\[\\neg\\mathrm\{stuck\}\],\\,n\_\{\\mathrm\{anti\}\},\\,\\mu,\\,\\alpha\)⊳\\trianglerightPhase 6

15:endif

16:endfor

17:return

\{xi:stucki\}\\\{x\_\{i\}:\\mathrm\{stuck\}\_\{i\}\\\}

### 3\.5GPU Execution Model

Every operation in Algorithm[1](https://arxiv.org/html/2606.26164#alg1)vectorizes over the sample indexii\. LetPPdenote the number of GPU parallel processors \(NVIDIA RTX 3080:P=8704P=8704CUDA cores\)\. For a population of sizeNN:

ForN=512N=512andd≤128d\\leq 128, all operations complete inO​\(1\)O\(1\)effective batches, yielding constant wall\-clock time per oscillation cycle regardless of population size up to GPU saturation\.

### 3\.6Hyperparameters

ChiSaohas several hyperparameters\. Table[1](https://arxiv.org/html/2606.26164#S3.T1)lists defaults and sensitivity range \(see Appendix[C](https://arxiv.org/html/2606.26164#A3)for full sensitivity analysis\)\.

Table 1:ChiSao hyperparameters with defaults\. All defaults were fixed before benchmarking and not tuned per\-problem\.

## 4Theoretical Analysis

### 4\.1Log\-Concave Case

###### Proposition 1\(Convergence for log\-concaveff\)\.

Letf:Θ→ℝf:\\Theta\\to\\mathbb\{R\}be strictly log\-concave on a convex compact domainΘ⊂ℝd\\Theta\\subset\\mathbb\{R\}^\{d\}with unique maximumθ∗∈int​\(Θ\)\\theta^\{\*\}\\in\\mathrm\{int\}\(\\Theta\)\. Then for any initial population\{xi\}i=1N\\\{x\_\{i\}\\\}\_\{i=1\}^\{N\}, after one oscillation cycle with sufficientnconvn\_\{\\mathrm\{conv\}\}, at least one sample satisfies‖xi−θ∗‖<ϵ\\\|x\_\{i\}\-\\theta^\{\*\}\\\|<\\epsilonwith probability 1, and that sample is marked stuck\.

###### Proof sketch\.

For strictly log\-concaveff, every stationary point offfonΘ\\Thetais the unique global maximumθ∗\\theta^\{\*\}\[[19](https://arxiv.org/html/2606.26164#bib.bib22)\]\. L\-BFGS with exact gradients converges to the unique stationary point for strongly convex objectives\[[12](https://arxiv.org/html/2606.26164#bib.bib17),[17](https://arxiv.org/html/2606.26164#bib.bib27)\]; the log\-concave case follows by the same argument applied to−f\-f\. The stick condition‖∇f​\(xi\)‖∞<ϵgrad\\\|\\nabla f\(x\_\{i\}\)\\\|\_\{\\infty\}<\\epsilon\_\{\\mathrm\{grad\}\}is satisfied at convergence for sufficiently small step sizes\. The quality thresholdf​\(xi\)≥fmax\+log⁡δf\(x\_\{i\}\)\\geq f\_\{\\max\}\+\\log\\deltais satisfied at the global maximum withfmax=f​\(θ∗\)f\_\{\\max\}=f\(\\theta^\{\*\}\)\. ∎

### 4\.2Multimodal Case: Coverage Analysis

For non\-log\-concaveffwithK∗K^\{\*\}significant modes\{θk∗\}k=1K∗\\\{\\theta^\{\*\}\_\{k\}\\\}\_\{k=1\}^\{K^\{\*\}\}, we analyze the probability that all modes are discovered\.

Letℬk=\{θ:θk∗=arg​maxθ′⁡f​\(θ′\),L​\-​BFGS​\(θ\)→θk∗\}\\mathcal\{B\}\_\{k\}=\\\{\\theta:\\theta^\{\*\}\_\{k\}=\\operatorname\*\{arg\\,max\}\_\{\\theta^\{\\prime\}\}f\(\\theta^\{\\prime\}\),\\mathrm\{L\\text\{\-\}BFGS\}\(\\theta\)\\to\\theta^\{\*\}\_\{k\}\\\}denote the*basin of attraction*of modekkunder L\-BFGS, with volumeVk=vol​\(ℬk∩Θ\)V\_\{k\}=\\mathrm\{vol\}\(\\mathcal\{B\}\_\{k\}\\cap\\Theta\)and total basin volumeVbasin=∑kVkV\_\{\\mathrm\{basin\}\}=\\sum\_\{k\}V\_\{k\}\.

###### Proposition 2\(Mode coverage probability\)\.

Let the initial population\{xi\}i=1N\\\{x\_\{i\}\\\}\_\{i=1\}^\{N\}be drawn i\.i\.d\. uniformly fromΘ\\Theta, and supposeVbasin/vol​\(Θ\)=ρ\>0V\_\{\\mathrm\{basin\}\}/\\mathrm\{vol\}\(\\Theta\)=\\rho\>0\(basins cover aρ\\rhofraction of the domain\)\. The probability that modekkis not discovered in the first convergence phase is:

P​\(miss mode​k\)=\(1−Vkvol​\(Θ\)\)N≤\(1−ρmin\)NP\(\\text\{miss mode \}k\)=\\left\(1\-\\frac\{V\_\{k\}\}\{\\mathrm\{vol\}\(\\Theta\)\}\\right\)^\{N\}\\leq\(1\-\\rho\_\{\\min\}\)^\{N\}\(7\)whereρmin=mink⁡Vk/vol​\(Θ\)\\rho\_\{\\min\}=\\min\_\{k\}V\_\{k\}/\\mathrm\{vol\}\(\\Theta\)\. ForN≥log⁡\(K∗/α\)/ρminN\\geq\\log\(K^\{\*\}/\\alpha\)/\\rho\_\{\\min\}, allK∗K^\{\*\}modes are discovered with probability at least1−α1\-\\alphaby union bound\.

### 4\.3GPU Complexity

###### Proposition 3\(Wall\-clock complexity\)\.

LetPPdenote GPU parallelism,NNpopulation size,dddimension,K∗K^\{\*\}number of modes, andTfT\_\{f\}the wall\-clock cost of a single function evaluation\. The total wall\-clock cost ofChiSaowithnoscn\_\{\\mathrm\{osc\}\}oscillations is:

TChiSao=O​\(nosc⋅\(nconv\+nanti\+ncloud​K\)⋅⌈N​dP⌉⋅Tf\)T\_\{\\textsc\{ChiSao\}\}=O\\\!\\left\(n\_\{\\mathrm\{osc\}\}\\cdot\\left\(n\_\{\\mathrm\{conv\}\}\+n\_\{\\mathrm\{anti\}\}\+n\_\{\\mathrm\{cloud\}\}K\\right\)\\cdot\\left\\lceil\\frac\{Nd\}\{P\}\\right\\rceil\\cdot T\_\{f\}\\right\)\(8\)In the GPU\-saturated regime \(N​d≤PNd\\leq P\), this reduces toO​\(nosc⋅nconv⋅Tf\)O\(n\_\{\\mathrm\{osc\}\}\\cdot n\_\{\\mathrm\{conv\}\}\\cdot T\_\{f\}\)—independent of bothNNanddd\.

Increasing the population sizeNNtherefore costs no additional wall\-clock time until the GPU is saturated \(N\>P/dN\>P/d\)\. Ford=64d=64on the RTX 3080 \(P≈8704P\\approx 8704\) the threshold isN≈136N\\approx 136; the defaultN=512N=512is at the threshold for this dimension\. Ford=2d=2, populations as large asN=4352N=4352are free\.

### 4\.4Limitations of the Analysis

The convergence guarantee \(Proposition[1](https://arxiv.org/html/2606.26164#Thmproposition1)\) requires log\-concavity, which fails for all genuinely multimodal functions\. The coverage analysis \(Proposition[2](https://arxiv.org/html/2606.26164#Thmproposition2)\) requires that basin volumes are non\-negligible and that the initial population is uniform; both assumptions can fail for pathological functions \(exponentially small basins, degenerate geometry\)\. The anti\-convergence momentum dynamics \(Phase 6\) are not analyzed; we have no guarantee that they increase coverage beyond the first\-pass uniform bound\. We claim no more\. The non\-log\-concave case rests on the empirical evidence of Section[5](https://arxiv.org/html/2606.26164#S5)\.

## 5Experiments

All experiments run on a single NVIDIA RTX 3080 Laptop GPU \(8 GB VRAM, 48 SMs, 8704 CUDA cores\)\. CPU baselines run on Intel Core i9\-12900H \(14 cores, 20 threads\)\.ChiSaois implemented in Python using CuPy\[[18](https://arxiv.org/html/2606.26164#bib.bib18)\]for GPU execution, with NumPy\[[7](https://arxiv.org/html/2606.26164#bib.bib20)\]fallback\. All hyperparameters use the defaults in Table[1](https://arxiv.org/html/2606.26164#S3.T1); no per\-problem tuning was performed\. Every benchmark function is evaluated value\-only:ChiSaocomputes its gradients by finite differences, so every wall\-clock time reported below is the derivative\-free worst case\[[22](https://arxiv.org/html/2606.26164#bib.bib28)\]\. Analytic or automatic\-differentiation gradients, when available, only reduce it\.

### 5\.1Benchmark Functions

We evaluate on all 42 standard test functions from the Simon Fraser University optimization benchmark suite\[[15](https://arxiv.org/html/2606.26164#bib.bib19)\]; seeJamil and Yang \[[8](https://arxiv.org/html/2606.26164#bib.bib26)\]for the wider catalogue andLiet al\.\[[11](https://arxiv.org/html/2606.26164#bib.bib29)\]for the CEC niching competition suite, spanning qualitatively distinct landscape types\. Functions are organized into four groups \(Table[2](https://arxiv.org/html/2606.26164#S5.T2)\): scalable multimodal \(Group A, tested acrossd∈\{2,4,8,16,32,64\}d\\in\\\{2,4,8,16,32,64\\\}\), scalable bowl and valley \(Group B, same dimensions\), fixed\-2D multimodal \(Group C\), and fixed low\-dimensional structured functions \(Group D\)\. We use 10 independent trials per condition with independent random initial populations of sizeN=200N=200\.

Table 2:All 42 benchmark functions, grouped by landscape character\. Group A and B functions are tested atd∈\{2,4,8,16,32,64\}d\\in\\\{2,4,8,16,32,64\\\}\. Group C functions atd=2d=2\. Group D functions at their canonical dimension\(s\)\.FunctionddtestedCharacterGroup A: Scalable multimodalRastrigin64Dense sinusoidal grid;≈\(2⋅5\.12/0\.5\)d\\approx\(2\\cdot 5\.12/0\.5\)^\{d\}local optimaAckley64Flat outer region; exponentially isolated central basinSchwefel64Deceptive: global optimum far from nearest competitorGriewank64Widespread local minima with multiplicative couplingLevy64Sinusoidally structured local minimaStyblinski\-Tang64Multiple local minima, quartic separableMichalewicz64d\!d\!local minima, steep channel ridges \(m=10m=10\)Group B: Scalable bowl / valleySphere64Strictly convex unimodalSum of Diff\. Powers64Unimodal, variable\-exponent separableRosenbrock64Narrow parabolic valleyZakharov64Unimodal, no spurious local minimaDixon\-Price64Analytic global min atxi=2−\(2i−2\)/2ix\_\{i\}=2^\{\-\(2^\{i\}\-2\)/2^\{i\}\}Trid64Quadratic bowl; analytic min atxi=i​\(D\+1−i\)x\_\{i\}=i\(D\+1\-i\), bounds\[−d2,d2\]\[\-d^\{2\},d^\{2\}\]Rotated Hyper\-Ellipsoid64Unimodal bowl with increasing ellipsoidal ridgeSum Squares64Unimodal weighted separableGroup C: Fixed\-2D multimodalEasom2Tiny isolated basin in large flat domainCross\-in\-Tray2Four equivalent global optimaDrop Wave2Wave\-like multimodal structureEggholder2Many sinusoidal local optima, difficult landscapeHolder Table2Four equivalent global optimaSchaffer N\.22Near\-circular oscillating gradientSchaffer N\.42Two global optima near axesLevy N\.132Fixed\-2D sinusoidal, asymmetricLangermann2Multimodal with damped oscillationsDe Jong N\.5225 shallow wellsShubert218 equivalent global optimaBukin N\.62Narrow ridge along parabolaBohachevsky2Multimodal with oscillating cosine termsGroup D: Fixed low\-ddstructured \(all methods 100%; see text\)Three\-hump Camel2Three local minima, one globalSix\-hump Camel2Six local minima, two globalBooth2Unimodal plateMatyas2Nearly flat unimodalMcCormick2Simple bimodalBeale2Saddle\-like, large flat regionsBranin2Three equivalent global optimaGoldstein\-Price2Dense multimodalHartmann 33Four local optimaHartmann 44Four local optimaHartmann 66Six local optimaShekel4Parameterised multimodal wellsColville4Four\-variable valleyPowell4–16Quadratic valley, multiple coupled axes
### 5\.2Metrics

We report two metrics\.Mode recovery rate: fraction of trials in which the global optimum is located within a function\-specificL∞L\_\{\\infty\}tolerance \(Table[2](https://arxiv.org/html/2606.26164#S5.T2)\)\.Mean wall\-clock time: total runtime in seconds averaged over 10 trials\.

### 5\.3Baselines

We compareChiSaoagainst three standard baselines\.Differential Evolution\(DE\): SciPy implementation\[[26](https://arxiv.org/html/2606.26164#bib.bib21)\]with default strategy\.Basin\-Hopping\(BH\): SciPy implementation with step\-size adaptation\.CMA\-ES:cmapackage with multi\-start restarts\.ChiSaois evaluated with two seeding strategies:random\(uniform random initialization\) andcarry\_tiger\(Carry Tiger to Mountain: structured ray\-based initialization from domain vertices, edges, and faces, matching the CarryTiger seeding of the SunBURST inference pipeline\[[29](https://arxiv.org/html/2606.26164#bib.bib23)\]\)\. All methods use matched function evaluation budgets\.

![Refer to caption](https://arxiv.org/html/2606.26164v1/fig2_success_rate_updated.png)Figure 1:Mean mode recovery rate vs\. dimensiondd, averaged across all Group A and B functions\.ChiSao\(both seeders, solid\) degrades gracefully; baselines \(dashed\) collapse sharply atd≥8d\\geq 8\. The mean is depressed for all methods by Schwefel, which no method solves aboved=8d=8\.
### 5\.4Mode Recovery Results

Tables[3](https://arxiv.org/html/2606.26164#S5.T3)–[5](https://arxiv.org/html/2606.26164#S5.T5)report mode recovery rates ford∈\{2,8,32,64\}d\\in\\\{2,8,32,64\\\}\(10 trials per condition\)\. Group D functions are omitted from these tables: all five methods achieve 100% recovery on every Group D function at its canonical dimension, with no discriminating information\. Full Group D results appear in the supplementary material\.

##### Scalable multimodal functions \(Group A\)\.

The dominant result is consistent 100% recovery by bothChiSaoseeders on Rastrigin, Levy, Styblinski\-Tang, and Michalewicz throughd=64d=64\. On Rastrigin, DE achieves 90% atd=2d=2but collapses entirely atd≥8d\\geq 8; BH and CMA\-ES follow the same pattern\. On Levy and Styblinski\-Tang, DE and CMA\-ES remain competitive throughd=8d=8before collapsing, while bothChiSaoseeders maintain perfect recovery\. Michalewicz—withd\!d\!local minima and steep ridges—is the only scalable multimodal function where all five methods achieve 100% at every dimension; the discriminator there is wall\-clock time \(Section[5\.5](https://arxiv.org/html/2606.26164#S5.SS5)\)\.

The Ackley result reveals a structural difference between the two seeding strategies\. Random seeding achieves 80% atd=4d=4and collapses entirely atd≥8d\\geq 8; carry\_tiger maintains 100% fromd=2d=2throughd=64d=64\. Ackley’s nearly\-flat outer region provides essentially no gradient signal toward the central basin for uniformly\-random samples, while the structured ray seeding cuts directly through the flat region\. Seeding is therefore not a tuning detail\. It is part of the algorithm, and the harder the landscape, the more it matters\. DE also reaches 100% throughd=8d=8before collapsing atd≥16d\\geq 16\.

On Griewank, carry\_tiger is the*only*method achieving 100% atd=2d=2andd=4d=4\(random seeding: 10% and 0% respectively; DE: 60% and 10%; CMA\-ES: 50% and 0%\)\. Fromd≥8d\\geq 8onward, bothChiSaoseeders fully recover, and DE partially recovers atd≥32d\\geq 32\(80%\)\. The carry\_tiger advantage at low\-ddGriewank shares the same mechanism as Ackley: product\-coupled local minima produce a gradient field that uniform random samples fail to navigate\.

Schwefel is the principal failure mode forChiSao\. Both seeders fail atd≥16d\\geq 16\(reporting no peaks, denoted†\\dagger\); carry\_tiger achieves only 10% atd=8d=8and random fails entirely\. The quality gate rejects all converged samples because Schwefel’s deceptive landscape drives convergence to secondary peaks whose log\-likelihood relative to the population maximum falls below the thresholdlog⁡0\.1\\log 0\.1\. DE achieves 80% atd=8d=8before also collapsing atd≥16d\\geq 16; no method succeeds beyondd=8d=8\. This failure is outsideChiSao’s design envelope, discussed in Section[6](https://arxiv.org/html/2606.26164#S6)\.

##### Scalable bowl and valley functions \(Group B\)\.

ChiSaoachieves 100% recovery on Sphere, Zakharov, Dixon\-Price, Rotated Hyper\-Ellipsoid, and Sum Squares throughd=64d=64\. CMA\-ES fails on Sphere atd≥16d\\geq 16\(despite it being strictly convex\) due to tolerance miscalibration in the multi\-restart configuration, and fails on Rotated Hyper\-Ellipsoid fromd≥16d\\geq 16for the same reason\.

On Rosenbrock, both seeders fail atd=2d=2\(20% and 50% respectively\) because the narrow parabolic valley floor creates a stationary manifold that fails the gradient\-norm stick condition, but both recover fully atd≥4d\\geq 4\. DE and BH succeed atd=2d=2; CMA\-ES collapses fromd=8d=8\.

Trid at high dimension is a quality\-gate failure: the broad quadratic bowl with domain\[−d2,d2\]\[\-d^\{2\},d^\{2\}\]provides insufficient gradient contrast atd≥32d\\geq 32\. Random seeding achieves 10% atd=32d=32and 0% atd=64d=64; carry\_tiger achieves 30% atd=32d=32and 0% atd=64d=64\. DE, BH, and CMA\-ES succeed at all dimensions since they apply no quality gate\.

Sum of Different Powers shows a subtler pattern:ChiSaorandom drops to 60% atd=64d=64while carry\_tiger remains at 100%, confirming that structured ray coverage provides a meaningful advantage at the largest dimensions\.

##### Fixed\-2D multimodal functions \(Group C\)\.

Table[5](https://arxiv.org/html/2606.26164#S5.T5)reports results for all 13 Group C functions\. BothChiSaoseeders achieve 100% on all but two\. Easom \(20% random, 60% CT\) presents an isolated basin covering roughly10−810^\{\-8\}of the domain area; carry\_tiger’s structured rays partially compensate but cannot guarantee coverage\. Bukin N\.6 \(0% both seeders\) fails entirely due to its narrow parabolic ridge, the same quality\-gate mechanism as Trid at highdd; DE, BH, and CMA\-ES all achieve 100%\. Drop Wave is the only function whereChiSaooutperforms a pure baseline: both seeders achieve 100% while DE reaches only 60% and CMA\-ES 40%\.

Table 3:Mode recovery rate \(%\) on Group A \(scalable multimodal\) functions,d∈\{2,8,32,64\}d\\in\\\{2,8,32,64\\\}, 10 trials\.Bold: best at each dimension\.†\\dagger: quality\-gate failure \(no peaks reported; see Section[6](https://arxiv.org/html/2606.26164#S6)\)\.Table 4:Mode recovery rate \(%\) on Group B \(scalable bowl and valley\) functions,d∈\{2,8,32,64\}d\\in\\\{2,8,32,64\\\}, 10 trials\.†\\dagger: quality\-gate failure on broad low\-contrast landscape \(see Section[6](https://arxiv.org/html/2606.26164#S6)\)\.Table 5:Mode recovery rate \(%\) on Group C \(fixed\-2D multimodal\) functions, 10 trials\.

### 5\.5Wall\-Clock Scaling

![Refer to caption](https://arxiv.org/html/2606.26164v1/fig3_wall_time_updated.png)Figure 2:Mean wall\-clock time vs\. dimensionddon a log scale, averaged over all functions where the method achieves\>0%\>0\\%recovery\. Fitted scaling exponents:ChiSao∝d0\.14\\propto d^\{0\.14\}–d0\.15d^\{0\.15\}\(GPU batch parallelism\); DE∝d0\.66\\propto d^\{0\.66\}; BH∝d0\.73\\propto d^\{0\.73\}; CMA\-ES∝d0\.53\\propto d^\{0\.53\}\.Table[6](https://arxiv.org/html/2606.26164#S5.T6)reports mean wall\-clock times for six representative functions covering the qualitatively distinct outcomes in the benchmark: functions where onlyChiSaosucceeds at highdd\(Rastrigin, Ackley, Levy, Styblinski\-Tang\), functions where all methods succeed \(Michalewicz\), and functions where all methods succeed but wall\-clock diverges dramatically \(Rotated Hyper\-Ellipsoid\)\. Times are shown for methods with\>0%\>0\\%recovery at that dimension; exact recovery rates are given in Tables[3](https://arxiv.org/html/2606.26164#S5.T3)–[4](https://arxiv.org/html/2606.26164#S5.T4)\.

The central performance claim is not accuracy alone but the combination of maintained accuracy and near\-constant wall\-clock across dimension\. On Rastrigin,ChiSao\(random\) runs in 1\.5–4\.2 s acrossd=2d=2tod=64d=64; all baselines require increasing time for 0% recovery fromd≥8d\\geq 8onward\. On Styblinski\-Tang,ChiSao\(random\) completes in 2\.8–5\.2 s across the full range while DE and CMA\-ES collapse atd≥32d\\geq 32\.

On Levy atd=64d=64,ChiSao\(random\) achieves 100% recovery in 25\.3 s; all baselines fail \(0%\) despite requiring 44\.5–315 s\. The combination of maintained recovery and faster absolute wall\-clock is the GPU dividend: work that scales with dimension on CPU is parallelized across the batch\.

On Michalewicz—the only scalable multimodal function where all methods achieve 100%—wall\-clock is the sole discriminator\. Atd=64d=64,ChiSao\(random\) completes in 23\.8 s and carry\_tiger in 29\.2 s; BH requires 805\.4 s \(∼34×\{\\sim\}34\\timesslower\), CMA\-ES 46\.5 s \(2×2\\times\), and DE 36\.1 s \(1\.5×1\.5\\times\)\. The34×34\\timesBH–ChiSaoratio atd=64d=64represents the GPU dividend on a problem where quality is not in question\.

Rotated Hyper\-Ellipsoid, a smooth unimodal function, provides the most extreme wall\-clock contrast: atd=64d=64, all four methods achieving 100% require 50\.2 s \(ChiSaorandom\), 99\.0 s \(CT\), 312\.5 s \(DE\), and 1946\.4 s \(BH\), a39×39\\timesgap betweenChiSaoand BH\. Since no multimodal complexity is involved, the entire speedup is attributable to GPU batch parallelism versus sequential CPU function evaluation\.

On Ackley, the carry\_tiger wall\-clock shows a non\-monotone pattern: 8\.6 s \(d=2d=2\), 7\.0 s \(d=8d=8\), 4\.6 s \(d=32d=32\), 23\.1 s \(d=64d=64\)\. The decrease fromd=2d=2tod=32d=32reflects early oscillation\-cycle termination when carry\_tiger seeding places samples close to the global basin, reducing the number of exploration phases required\. Thed=64d=64increase reflects additional L\-BFGS iterations needed to converge in the high\-dimensional flat landscape\.

Table 6:Mean wall\-clock time \(seconds\) over 10 trials atd∈\{2,8,32,64\}d\\in\\\{2,8,32,64\\\}\. Entries marked “—” indicate 0% recovery at that dimension\. Times for methods with partial recovery \(<100%<100\\%\) reflect only the trials that ran to completion; exact rates appear in Tables[3](https://arxiv.org/html/2606.26164#S5.T3)–[4](https://arxiv.org/html/2606.26164#S5.T4)\.ChiSao\(CT\) denotes carry\_tiger seeder\.
### 5\.6Ablation Study

![Refer to caption](https://arxiv.org/html/2606.26164v1/fig1_ablation.png)Figure 3:Ablation heatmap: mode recovery rate atd=8d=8for each phase removal, under random \(left\) and carry\_tiger \(right\) seeding\. Green = 100%, red = 0%\. Removing deduplication \(Ph\. 3\) is universally catastrophic\. Removing anti\-convergence \(Ph\. 6\) improves Schwefel atd=8d=8under CT seeding, but this effect does not persist at higher dimensions \(see text\)\.Table 7:Ablation study: mode recovery rate \(%\) on Group A functions atd=8d=8, 10 trials\. Each row disables one structural component\. Two seeding strategies shown: random \(left\) and carry\_tiger \(right, italicised\)\.Bold: degraded relative to fullChiSaowith the same seeder\.Three results dominate the ablation \(Table[7](https://arxiv.org/html/2606.26164#S5.T7)\), with thed=8d=8findings extended by additional ablations atd∈\{16,32,64\}d\\in\\\{16,32,64\\\}under carry\_tiger seeding\.

##### Deduplication is not optional\.

Disabling Phase 3 \(no deduplication\) causes 0% recovery on all seven functions under both seeders at every tested dimension\. Without deduplication, the reseeding phase \(Phase 4\) replenishes consumed slots from known peaks without removing them from the stuck set; the population rapidly saturates at the first discovered mode and exploration terminates\. Deduplication is the only component whose removal is universally and persistently catastrophic\.

##### FullChiSaois the most robust configuration across dimensions\.

Across all tested dimensions \(d∈\{8,16,32,64\}d\\in\\\{8,16,32,64\\\}\), fullChiSaowith carry\_tiger seeding achieves a stable 6/7 success rate, recovering all Group A functions except Schwefel consistently\. No other configuration maintains this stability\. Plain multistart \(all off\) degrades to 5/7 atd≥32d\\geq 32as Levy recovery collapses\. Schwefel remains the persistent exception at all dimensions and all configurations: its deceptive landscape places the global basin far from secondary peaks, and the quality\-gate interacts unfavourably with the gradient structure at higher dimensions regardless of phase configuration\.

##### Anti\-convergence effects are dimension\-dependent\.

Atd=8d=8under carry\_tiger, disabling Phase 6 raises Schwefel recovery from 0% to 100%, appearing to make it the optimal minimal configuration\. However, this advantage is dimension\-specific: atd=16d=16, Levy recovery under the Ph\. 6\-off configuration drops to 20%; atd=32d=32, Levy fails entirely while Schwefel recovers only 60%; atd=64d=64, both Levy and Schwefel fail \(0%\) under Ph\. 6\-off\. Anti\-convergence is what sustains Levy recovery at scale\. The mechanism is as follows: on Schwefel atd=8d=8, the Hands Like Clouds phase navigates samples toward the global region but the anti\-convergence momentum then displaces them into secondary peaks that trigger quality\-gate failure; disabling anti\-convergence fixes this at low dimension\. At higher dimensions, however, basin volumes contract and anti\-convergence becomes necessary for Levy to direct samples away from shallow traps during later oscillation passes\. Thed=8d=8ablation result is thus not a stable operating point—it resolves one failure mode while introducing another that emerges at scale\. Under random seeding, anti\-convergence is additionally critical for Rastrigin at all dimensions \(0% without it\), reinforcing that it remains a necessary component of the full algorithm\.

### 5\.7Noise Robustness

We testChiSao’s mode detection reliability whenffis corrupted by additive Gaussian noise:f~​\(θ\)=f​\(θ\)\+ε\\tilde\{f\}\(\\theta\)=f\(\\theta\)\+\\varepsilon,ε∼𝒩​\(0,σnoise2\)\\varepsilon\\sim\\mathcal\{N\}\(0,\\sigma\_\{\\mathrm\{noise\}\}^\{2\}\)\. We use a 2\-mode Gaussian mixture ind=6d=6with peaks at±2​e^1\\pm 2\\hat\{e\}\_\{1\}, mode separation4​σ4\\sigma\.

These results are available from the SunBURST benchmark suite\[[29](https://arxiv.org/html/2606.26164#bib.bib23)\]; we reproduce them here for completeness\.

Table 8:Mode detection reliability under likelihood noise \(d=6d=6, 2 peaks at±2​e^1\\pm 2\\hat\{e\}\_\{1\}, 10 trials per level\)\. Both peaks are correctly identified in all 70 trials across all noise levels\.The key result is 100% mode detection even atσnoise=1\.0\\sigma\_\{\\mathrm\{noise\}\}=1\.0, where noise amplitude equals signal scale\. Peak\-location error degrades proportionally to noise as expected from Hessian estimation uncertainty, but mode detection is never compromised\. Zero spurious detections throughσ=0\.5\\sigma=0\.5\. This robustness follows fromChiSao’s gradient\-based stick detection: the convergence\-anticonvergence cycle is driven by sign\-consistent gradient estimates averaged over the trajectory rather than evaluated at a single noisy point\.

## 6Discussion

The noise result is the one that matters \(Table[8](https://arxiv.org/html/2606.26164#S5.T8)\)\. Mode detection holds at 100% even when the noise amplitude equals the signal scale\. This is not luck: trajectory\-averaged gradients absorb the per\-evaluation noise that point\-evaluation methods cannot tolerate\. The oscillation cycle, the freeze\-and\-explore asymmetry, and the reseeding strategies do the rest\.

##### Relation to SunBURST\.

ChiSaowas developed as the mode\-discovery engine for SunBURST\[[29](https://arxiv.org/html/2606.26164#bib.bib23)\], a GPU\-accelerated Bayesian evidence calculator\. In that context, ChiSao receives initial populations from a ray\-casting seeding strategy \(CarryTiger\) and its output modes are passed to a Laplace\-approximation evidence integrator \(BendTheBow\)\. The present paper establishes ChiSao as a standalone contribution independent of that pipeline, applicable to any black\-box optimization problem with gradient access\.

##### Limitations\.

ChiSaoassumes a smooth objective: it estimates gradients by finite differences \(the default, used for every result here\) or accepts them analytically, so genuinely non\-differentiable or discrete problems are out of scope\. The Hands Like Clouds smoothing scaleσ\\sigmais auto\-estimated but may require manual tuning on strongly heterogeneous landscapes\. The anti\-convergence momentum step sizeα\\alphais sensitive on functions with very small or very large gradient magnitudes; an adaptive line\-search would improve robustness\. The theoretical analysis \(Section[4](https://arxiv.org/html/2606.26164#S4)\) covers only the log\-concave and first\-oscillation\-pass cases; a full convergence theory for the oscillation dynamics remains open\.

##### Failure modes\.

The quality\-gate threshold \(δ=0\.1\\delta=0\.1, i\.e\., modes within one log\-decade of the global maximum\) is the primary source ofChiSao’s failures on this benchmark\. On Schwefel, the deceptive landscape drives convergence to secondary peaks whose likelihood falls below threshold; on Trid at highddand Bukin N\.6, the broad low\-contrast bowl and narrow ridge, respectively, produce converged samples that pass neither the gradient\-norm nor the likelihood test\. In all three cases, the failure is a property of the quality gate’s interaction with the landscape, not of the search dynamics\. Baselines without quality gates \(DE, BH, CMA\-ES\) succeed whereChiSaofails on these functions, at the cost of reporting spurious peaks on harder multimodal problems\.ChiSaounderperforms additionally on functions with exponentially many equal\-height modes: the deduplication and reseeding logic does not scale toO​\(2d\)O\(2^\{d\}\)distinct modes\.

##### Future work\.

Since all reported results already use finite\-difference gradients, analytic or automatic\-differentiation gradients are a direct speedup lever: they remove the2​d​N2dNfinite\-difference evaluations per L\-BFGS step, though the GPU already absorbs them inO​\(1\)O\(1\)wall\-clock while2​d​N<P2dN<P\. On the theoretical side the missing piece is a Markov\-chain analysis of basin\-crossing probability under the anti\-convergence momentum dynamics, which would close the gap between Proposition[1](https://arxiv.org/html/2606.26164#Thmproposition1)and the empirical multi\-modal results\.

## Data Availability

## Code Availability

ChiSaois released as a standalone open\-source Python package \(MIT License\) on PyPI, installable viapip install chisao\. Source code, the test suite, and the full SFU benchmark suite are available at[https://github\.com/beastraban/chisao](https://github.com/beastraban/chisao)\.

## Competing Interests

The author declares no competing interests\.

## Appendix ATheL∞L\_\{\\infty\}Metric for Deduplication

TheL∞L\_\{\\infty\}\(Chebyshev\) distance‖x−y‖∞=maxi⁡\|xi−yi\|\\\|x\-y\\\|\_\{\\infty\}=\\max\_\{i\}\|x\_\{i\}\-y\_\{i\}\|has two advantages overL2L\_\{2\}for high\-dimensional deduplication\. First, it is dimension\-independent in expectation: for uniform random points in\[0,1\]d\[0,1\]^\{d\},𝔼​\[‖x−y‖∞\]≈d/\(d\+1\)→1\\mathbb\{E\}\[\\\|x\-y\\\|\_\{\\infty\}\]\\approx d/\(d\+1\)\\to 1, whileL2L\_\{2\}distance diverges asd/6\\sqrt\{d/6\}\. A fixed thresholdϵdup\\epsilon\_\{\\mathrm\{dup\}\}inL∞L\_\{\\infty\}therefore has consistent semantics across dimensions\. Second, theL∞L\_\{\\infty\}ball\{y:‖x−y‖∞≤r\}\\\{y:\\\|x\-y\\\|\_\{\\infty\}\\leq r\\\}is a hypercube of side2​r2r, enabling efficient spatial hashing forO​\(N\)O\(N\)deduplication rather thanO​\(N2\)O\(N^\{2\}\)pairwise comparison\.

## Appendix BHands Like Clouds: Smoothed Gradient Analysis

The smoothed functionfσ​\(θ\)=𝔼z∼𝒩​\(0,Id\)​\[f​\(θ\+σ​z\)\]f\_\{\\sigma\}\(\\theta\)=\\mathbb\{E\}\_\{z\\sim\\mathcal\{N\}\(0,I\_\{d\}\)\}\[f\(\\theta\+\\sigma z\)\]satisfies:

∇fσ​\(θ\)=𝔼z∼𝒩​\(0,Id\)​\[∇f​\(θ\+σ​z\)\]\\nabla f\_\{\\sigma\}\(\\theta\)=\\mathbb\{E\}\_\{z\\sim\\mathcal\{N\}\(0,I\_\{d\}\)\}\[\\nabla f\(\\theta\+\\sigma z\)\]\(9\)by Leibniz’s rule \(assumingffis differentiable a\.e\.\)\. The Monte Carlo estimator \([5](https://arxiv.org/html/2606.26164#S3.E5)\) is unbiased\. For functions with multi\-scale structuref=fglobal\+flocalf=f\_\{\\mathrm\{global\}\}\+f\_\{\\mathrm\{local\}\}whereflocalf\_\{\\mathrm\{local\}\}has characteristic scaleλ<σ\\lambda<\\sigma, the smoothed gradient∇fσ\\nabla f\_\{\\sigma\}eliminatesflocalf\_\{\\mathrm\{local\}\}contributions \(their integral over𝒩​\(0,σ2​I\)\\mathcal\{N\}\(0,\\sigma^\{2\}I\)vanishes by symmetry for zero\-mean local oscillations\), revealing the global basin structure\.

Nesterov and Spokoiny \[[16](https://arxiv.org/html/2606.26164#bib.bib16)\]showed that for convexff, gradient descent onfσf\_\{\\sigma\}converges to withinO​\(σ2​L\)O\(\\sigma^\{2\}L\)of the optimum offf, whereLLis the Lipschitz constant of∇f\\nabla f\. For non\-convexff, the smoothing provides no convergence guarantee but empirically helps samples navigate past local barriers whose scale is smaller thanσ\\sigma\.

## Appendix CHyperparameter Sensitivity

Each row of Table[9](https://arxiv.org/html/2606.26164#A3.T9)varies one hyperparameter from its default while holding the rest fixed at the values of Table[1](https://arxiv.org/html/2606.26164#S3.T1)\. “Robust set” is the mean mode\-recovery rate across Rastrigin, Ackley, Levy, Griewank, Styblinski\-Tang, and Michalewicz,d=32d\{=\}32, averaged over both seeders \(random, carry\_tiger\), 30 trials per cell\. Schwefel is reported separately atd=8d\{=\}8because it sits at a noise floor for every configuration, as discussed below\.

Table 9:Hyperparameter sensitivity sweep\. Mode recovery rate \(%\), 30 trials per cell\.Three observations follow\. \(i\) On the six landscapes thatChiSaosolves at 100% with the defaults, recovery is invariant to every perturbation in the sweep with one exception\. \(ii\) That exception isnanti=10n\_\{\\mathrm\{anti\}\}=10: Levy collapses from 100% to 0% under both seeders\. Doubling the anti\-convergence step count over\-disperses samples past Levy’s narrow global basin during the escape phase, and the next L\-BFGS pass returns from outside the basin of attraction\. The defaultnanti=5n\_\{\\mathrm\{anti\}\}=5is on the safe side of this ceiling; this is the single non\-trivial constraint that the sweep imposes on the defaults\. \(iii\) Schwefel sits at a 5–13% recovery floor across every cell, including the default, with no swept hyperparameter raising it out of the noise band\. This confirms the failure\-mode diagnosis in Section[6](https://arxiv.org/html/2606.26164#S6): Schwefel’s failure is a property of the quality\-gate thresholdδ=0\.1\\delta=0\.1interacting with the deceptive landscape, not of the search\-dynamics parameters that the sweep covers\. A sensitivity sweep overδ\\deltawould close this loop, butδ\\deltais not currently exposed as a hyperparameter in thesticky\_handsinterface\.

A full per\-function breakdown of the sweep \(84 cells×\\times2 seeders\) is provided in the supplementary results file\.

## Appendix DBenchmark Function Definitions

For reproducibility, we give exact definitions used in experiments\. All functions are negated for maximization\.

##### Rastrigin\.

f​\(θ\)=−\[A​d\+∑i=1d\(θi2−A​cos⁡\(2​π​θi\)\)\]f\(\\theta\)=\-\\left\[Ad\+\\sum\_\{i=1\}^\{d\}\(\\theta\_\{i\}^\{2\}\-A\\cos\(2\\pi\\theta\_\{i\}\)\)\\right\],A=10A=10,θ∈\[−5\.12,5\.12\]d\\theta\\in\[\-5\.12,5\.12\]^\{d\}\. Global maximum atθ=0\\theta=0,f​\(0\)=0f\(0\)=0\.

##### Ackley\.

f​\(θ\)=−\[−a​exp⁡\(−b​d−1​∑θi2\)−exp⁡\(d−1​∑cos⁡\(c​θi\)\)\+a\+e\]f\(\\theta\)=\-\\left\[\-a\\exp\(\-b\\sqrt\{d^\{\-1\}\\sum\\theta\_\{i\}^\{2\}\}\)\-\\exp\(d^\{\-1\}\\sum\\cos\(c\\theta\_\{i\}\)\)\+a\+e\\right\],a=20a=20,b=0\.2b=0\.2,c=2​πc=2\\pi,θ∈\[−32\.768,32\.768\]d\\theta\\in\[\-32\.768,32\.768\]^\{d\}\.

##### Schwefel\.

f​\(θ\)=−\[418\.9829​d−∑i=1dθi​sin⁡\(\|θi\|\)\]f\(\\theta\)=\-\\left\[418\.9829d\-\\sum\_\{i=1\}^\{d\}\\theta\_\{i\}\\sin\(\\sqrt\{\|\\theta\_\{i\}\|\}\)\\right\],θ∈\[−500,500\]d\\theta\\in\[\-500,500\]^\{d\}\. Global maximum nearθi=420\.97\\theta\_\{i\}=420\.97\.

##### Styblinski\-Tang\.

f​\(θ\)=−12​∑i=1d\(θi4−16​θi2\+5​θi\)f\(\\theta\)=\-\\frac\{1\}\{2\}\\sum\_\{i=1\}^\{d\}\(\\theta\_\{i\}^\{4\}\-16\\theta\_\{i\}^\{2\}\+5\\theta\_\{i\}\),θ∈\[−5,5\]d\\theta\\in\[\-5,5\]^\{d\}\.

##### Michalewicz\.

f​\(θ\)=∑i=1dsin⁡\(θi\)​sin2​m⁡\(i​θi2/π\)f\(\\theta\)=\\sum\_\{i=1\}^\{d\}\\sin\(\\theta\_\{i\}\)\\sin^\{2m\}\(i\\theta\_\{i\}^\{2\}/\\pi\),m=10m=10,θ∈\[0,π\]d\\theta\\in\[0,\\pi\]^\{d\}\.

##### Trid\.

f​\(θ\)=−\[∑i=1d\(θi−1\)2−∑i=2dθi​θi−1\]f\(\\theta\)=\-\\left\[\\sum\_\{i=1\}^\{d\}\(\\theta\_\{i\}\-1\)^\{2\}\-\\sum\_\{i=2\}^\{d\}\\theta\_\{i\}\\theta\_\{i\-1\}\\right\],θ∈\[−d2,d2\]d\\theta\\in\[\-d^\{2\},d^\{2\}\]^\{d\}\. Global maximum atθi∗=i​\(d\+1−i\)\\theta\_\{i\}^\{\*\}=i\(d\+1\-i\)withf∗=d​\(d\+4\)​\(d−1\)/6f^\{\*\}=d\(d\+4\)\(d\-1\)/6\. Unimodal but increasingly ill\-conditioned withdd: Hessian eigenvalues span many orders of magnitude and the optimum coordinatesθi∗\\theta\_\{i\}^\{\*\}range fromddtod​\(d\+1\)2/4d\(d\+1\)^\{2\}/4\.

##### Bukin N\.6\.

f​\(x,y\)=−\[100​\|y−0\.01​x2\|\+0\.01​\|x\+10\|\]f\(x,y\)=\-\\left\[100\\sqrt\{\|y\-0\.01\\,x^\{2\}\|\}\+0\.01\\,\|x\+10\|\\right\],x∈\[−15,−5\]x\\in\[\-15,\-5\],y∈\[−3,3\]y\\in\[\-3,3\]\. Global maximum at\(−10,1\)\(\-10,1\)withf∗=0f^\{\*\}=0\. The minimum lies on the parabolic ridgey=0\.01​x2y=0\.01\\,x^\{2\}; the function is non\-differentiable along this ridge, and the\|⋅\|\\sqrt\{\|\\cdot\|\}gradient diverges as the ridge is approached\.

##### Gaussian mixture\.

f​\(θ\)=log​∑k=1Kwk​𝒩​\(θ;μk,σ2​I\)f\(\\theta\)=\\log\\sum\_\{k=1\}^\{K\}w\_\{k\}\\mathcal\{N\}\(\\theta;\\mu\_\{k\},\\sigma^\{2\}I\)withwk=1/Kw\_\{k\}=1/K,μk\\mu\_\{k\}random inΘ\\Theta,σ=1\\sigma=1,Θ=\[−10,10\]d\\Theta=\[\-10,10\]^\{d\}\.

## References

- \[1\]\(2020\)JAXNS: a high\-performance nested sampling package based on JAX\.arXiv preprint arXiv:2012\.15286\.Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px6.p1.1)\.
- \[2\]A\. Blake and A\. Zisserman\(1987\)Visual reconstruction\.MIT Press,Cambridge, MA\.Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px7.p1.1)\.
- \[3\]R\. Brits, A\. P\. Engelbrecht, and F\. van den Bergh\(2002\)A niching particle swarm optimizer\.Proceedings of the 4th Asia\-Pacific Conference on Simulated Evolution and Learning,pp\. 692–696\.Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px3.p1.1)\.
- \[4\]I\. Goodfellow, Y\. Bengio, and A\. Courville\(2016\)Deep learning\.MIT Press\.Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px6.p1.1)\.
- \[5\]B\. Hajek\(1988\)Cooling schedules for optimal annealing\.Mathematics of Operations Research13\(2\),pp\. 311–329\.Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px2.p1.1)\.
- \[6\]N\. Hansen and A\. Ostermeier\(2001\)Completely derandomized self\-adaptation in evolution strategies\.Evolutionary Computation9\(2\),pp\. 159–195\.Cited by:[§1](https://arxiv.org/html/2606.26164#S1.p2.1),[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px3.p1.1)\.
- \[7\]C\. R\. Harris, K\. J\. Millman, S\. J\. van der Walt,et al\.\(2020\)Array programming with NumPy\.Nature585\(7825\),pp\. 357–362\.Cited by:[§5](https://arxiv.org/html/2606.26164#S5.p1.1)\.
- \[8\]M\. Jamil and X\. Yang\(2013\)A literature survey of benchmark functions for global optimisation problems\.International Journal of Mathematical Modelling and Numerical Optimisation4\(2\),pp\. 150–194\.External Links:[Document](https://dx.doi.org/10.1504/IJMMNO.2013.055204),1308\.4008Cited by:[§5\.1](https://arxiv.org/html/2606.26164#S5.SS1.p1.2)\.
- \[9\]J\. Kennedy and R\. Eberhart\(1995\)Particle swarm optimization\.InProceedings of the IEEE International Conference on Neural Networks,Vol\.4,pp\. 1942–1948\.Cited by:[§1](https://arxiv.org/html/2606.26164#S1.p2.1),[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px3.p1.1)\.
- \[10\]S\. Kirkpatrick, C\. D\. Gelatt, and M\. P\. Vecchi\(1983\)Optimization by simulated annealing\.Science220\(4598\),pp\. 671–680\.Cited by:[§1](https://arxiv.org/html/2606.26164#S1.p2.1),[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px2.p1.1)\.
- \[11\]X\. Li, A\. Engelbrecht, and M\. G\. Epitropakis\(2013\)Benchmark functions for CEC’2013 special session and competition on niching methods for multimodal function optimization\.Technical reportEvolutionary Computation and Machine Learning Group, RMIT University\.Cited by:[§5\.1](https://arxiv.org/html/2606.26164#S5.SS1.p1.2)\.
- \[12\]D\. C\. Liu and J\. Nocedal\(1989\)On the limited memory BFGS method for large scale optimization\.Mathematical Programming45\(1\-3\),pp\. 503–528\.Cited by:[§3\.2\.1](https://arxiv.org/html/2606.26164#S3.SS2.SSS1.p1.3),[§4\.1](https://arxiv.org/html/2606.26164#S4.SS1.1.p1.8)\.
- \[13\]R\. Martí\(2003\)Multi\-start methods\.InHandbook of Metaheuristics,pp\. 355–368\.Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px4.p1.1)\.
- \[14\]H\. Mobahi and J\. W\. Fisher\(2015\)On the link between Gaussian homotopy continuation and convex envelopes\.InEnergy Minimization Methods in Computer Vision and Pattern Recognition \(EMMCVPR\),Lecture Notes in Computer Science, Vol\.8932,pp\. 43–56\.External Links:[Document](https://dx.doi.org/10.1007/978-3-319-14612-6%5F4)Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px7.p1.1)\.
- \[15\]M\. Molga and C\. Smutnicki\(2005\)Test functions for optimization needs\.Technical reportWrocław University of Technology\.Note:Available at[http://www\.zsd\.ict\.pwr\.wroc\.pl/files/docs/functions\.pdf](http://www.zsd.ict.pwr.wroc.pl/files/docs/functions.pdf)Cited by:[§5\.1](https://arxiv.org/html/2606.26164#S5.SS1.p1.2)\.
- \[16\]Y\. Nesterov and V\. Spokoiny\(2017\)Random gradient\-free minimization of convex functions\.Foundations of Computational Mathematics17\(2\),pp\. 527–566\.Cited by:[Appendix B](https://arxiv.org/html/2606.26164#A2.p2.8),[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px7.p1.1)\.
- \[17\]J\. Nocedal and S\. J\. Wright\(2006\)Numerical optimization\.2nd edition,Springer Series in Operations Research and Financial Engineering,Springer,New York\.External Links:ISBN 978\-0\-387\-30303\-1Cited by:[§3\.2\.1](https://arxiv.org/html/2606.26164#S3.SS2.SSS1.p1.3),[§4\.1](https://arxiv.org/html/2606.26164#S4.SS1.1.p1.8)\.
- \[18\]R\. Okuta, Y\. Unno, D\. Nishino, S\. Hido, and C\. Loomis\(2017\)CuPy: a NumPy\-compatible library for NVIDIA GPU calculations\.InProceedings of Workshop on Machine Learning Systems \(LearningSys\) in The Thirty\-first Annual Conference on Neural Information Processing Systems \(NeurIPS\),Cited by:[§5](https://arxiv.org/html/2606.26164#S5.p1.1)\.
- \[19\]A\. Prékopa\(1973\)On logarithmic concave measures and functions\.Acta Scientiarum Mathematicarum34,pp\. 335–343\.Cited by:[§4\.1](https://arxiv.org/html/2606.26164#S4.SS1.1.p1.8)\.
- \[20\]M\. Preuß\(2015\)Multimodal optimization by means of evolutionary algorithms\.Springer\.Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px3.p1.1)\.
- \[21\]D\. J\. Ram, T\. Sreenivas, and K\. G\. Subramaniam\(1996\)Parallel simulated annealing algorithms\.Journal of Parallel and Distributed Computing37\(2\),pp\. 207–212\.Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px2.p1.1)\.
- \[22\]L\. M\. Rios and N\. V\. Sahinidis\(2013\)Derivative\-free optimization: a review of algorithms and comparison of software implementations\.Journal of Global Optimization56\(3\),pp\. 1247–1293\.External Links:[Document](https://dx.doi.org/10.1007/s10898-012-9951-y)Cited by:[§5](https://arxiv.org/html/2606.26164#S5.p1.1)\.
- \[23\]T\. Salimans, J\. Ho, X\. Chen, S\. Sidor, and I\. Sutskever\(2017\)Evolution strategies as a scalable alternative to reinforcement learning\.arXiv preprint arXiv:1703\.03864\.External Links:1703\.03864Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px6.p1.1)\.
- \[24\]B\. Shahriari, K\. Swersky, Z\. Wang, R\. P\. Adams, and N\. de Freitas\(2016\)Taking the human out of the loop: a review of Bayesian optimization\.Proceedings of the IEEE104\(1\),pp\. 148–175\.External Links:[Document](https://dx.doi.org/10.1109/JPROC.2015.2494218)Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px5.p1.1)\.
- \[25\]R\. Storn and K\. Price\(1997\)Differential evolution – a simple and efficient heuristic for global optimization over continuous spaces\.Journal of Global Optimization11\(4\),pp\. 341–359\.Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px3.p1.1)\.
- \[26\]P\. Virtanen, R\. Gommers, T\. E\. Oliphant,et al\.\(2020\)SciPy 1\.0: fundamental algorithms for scientific computing in Python\.Nature Methods17,pp\. 261–272\.Cited by:[§5\.3](https://arxiv.org/html/2606.26164#S5.SS3.p1.1)\.
- \[27\]D\. J\. Wales and J\. P\. Doye\(1997\)Global optimization by basin\-hopping and the lowest energy structures of Lennard\-Jones clusters containing up to 110 atoms\.The Journal of Physical Chemistry A101\(28\),pp\. 5111–5116\.Cited by:[§1](https://arxiv.org/html/2606.26164#S1.p2.1),[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px1.p1.1),[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px8.p1.1)\.
- \[28\]D\. J\. Wales\(2003\)Energy landscapes\.Cambridge University Press\.Cited by:[§2](https://arxiv.org/html/2606.26164#S2.SS0.SSS0.Px1.p1.1)\.
- \[29\]I\. Wolfson\(2026\)SunBURST: deterministic GPU\-accelerated Bayesian evidence via mode\-centric Laplace integration\.arXiv preprint arXiv:2601\.19957\.External Links:2601\.19957Cited by:[§5\.3](https://arxiv.org/html/2606.26164#S5.SS3.p1.1),[§5\.7](https://arxiv.org/html/2606.26164#S5.SS7.p2.1),[§6](https://arxiv.org/html/2606.26164#S6.SS0.SSS0.Px1.p1.1)\.