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

1 Introduction↩︎

Multimodal optimization—finding all significant maxima of a function \(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 [1] and simulated annealing [2] 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 [3] and particle swarm optimization [4] 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 with \(P > 10^4\) parallel cores, a batch of \(N \leq P\) gradient 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.

ChiSao¹ 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; ChiSao advances 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 \(\nabla f\). Smoothing at scale \(\sigma\) erases sub-\(\sigma\) texture 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.

1.0.0.1 Contributions.

The contribution is fourfold: the oscillation cycle as a strategy for multimodal optimization with a formalized freeze-and-explore asymmetry (Section 3); convergence guarantees for log-concave targets and a coverage analysis for the multimodal case (Section 4); benchmarks against Differential Evolution, Basin-Hopping, and CMA-ES on all 42 SFU functions across \(d \in \{2,4,8,16,32,64\}\) (Section 5); and noise-robustness results showing 100% mode detection up to \(\sigma_{\mathrm{noise}} = 1.0\), the signal scale itself (Section 5.7).

2 Related Work↩︎

2.0.0.1 Basin-hopping.

Wales and Doye [1] 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 [5]. 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.

2.0.0.2 Simulated annealing.

Kirkpatrick et al.[2] introduced temperature-based stochastic acceptance, which provides theoretical convergence guarantees under slow cooling [6] but requires an exponentially slow schedule to avoid premature convergence in practice. Parallel implementations [7] distribute independent chains but do not share information, losing the benefit of population diversity.

2.0.0.3 Evolutionary and population methods.

CMA-ES [3] adapts a full covariance matrix to the population geometry, achieving state-of-the-art performance on unimodal and mildly multimodal problems. Niching extensions [8] attempt to maintain diversity but add significant complexity. Particle swarm optimization [4] uses social attraction toward the population best, which concentrates the population at a single basin in multimodal settings unless modified with repulsion terms [9]. Differential evolution [10] generates candidates by combining population members but shares CMA-ES’s tendency toward mode collapse.

2.0.0.4 Multistart methods.

Running multiple independent gradient-based optimizers from random initial points is the most common approach in practice [11]. 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.

2.0.0.5 Bayesian optimization.

Gaussian-process surrogate methods [12] 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.

2.0.0.6 GPU-parallel optimization.

GPU acceleration of optimization has focused primarily on stochastic gradient descent for deep learning [13], where the objective is unimodal and the bottleneck is data throughput rather than exploration. For black-box multimodal optimization, GPU acceleration is largely unexplored, with [14] the major exception—evolution strategies massively parallelized for reinforcement learning, but with a single-objective design that does not catalogue modes. [15] parallelized nested sampling on GPU (JAXNS) but the sequential compression structure of nested sampling limits GPU utilization to within-iteration batching. ChiSao appears to be the first black-box multimodal optimizer with full GPU utilization across the search strategy itself, not only inside individual iterations.

2.0.0.7 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” [16], [17]. [18] analyzed Gaussian smoothing for convex optimization. The Hands Like Clouds phase of ChiSao applies smoothing as one phase within a structured oscillation cycle, rather than as a global preprocessing step or annealing schedule.

2.0.0.8 Position relative to prior work.

ChiSao differs 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 [1] is the closest neighbour; it differs on all four.

3 Algorithm↩︎

3.1 Problem Statement↩︎

Let \(f: \Theta \to \mathbb{R}\) be a black-box function on a bounded domain \(\Theta \subset \mathbb{R}^d\), accessible through evaluations \(f(\theta)\); the gradients \(\nabla f(\theta)\) the algorithm uses are computed by finite differences, or supplied analytically when available. We seek the set of all significant modes: \[\mathcal{M}^* = \left\{ \theta^* \in \Theta : \nabla f(\theta^*) = 0,\; \nabla^2 f(\theta^*) \prec 0,\; f(\theta^*) \geq f_{\max} + \log \delta \right\}\] where \(f_{\max} = \max_\theta f(\theta)\) and \(\delta \in (0,1)\) is a quality threshold (default \(\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}^*\).

ChiSao takes as input a population \(\{x_i\}_{i=1}^N \subset \Theta\) of initial candidate points (arbitrary; provided by the caller), and returns the estimated set \(\hat{\mathcal{M}}^* \subseteq \{x_i\}\).

3.2 The Oscillation Cycle↩︎

ChiSao runs for \(n_{\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.

3.2.1 Phase 1: Convergence↩︎

All stuck masks are released (every sample is temporarily unfrozen) and the entire population undergoes \(n_{\mathrm{conv}}\) steps of batched L-BFGS [19], [20] toward local maxima of \(f\). All \(N\) samples optimize simultaneously in a single GPU batch: \[x_i \gets \mathrm{L\text{-}BFGS}(x_i, n_{\mathrm{conv}}) \quad \forall i = 1, \ldots, N \quad \text{(GPU-parallel)}\] The L-BFGS memory parameter is \(m = 10\) by default. Iteration count adapts to dimension: \(n_{\mathrm{conv}} = \max(10, 3\log_2 d)\).

3.2.2 Phase 2: Stick Detection↩︎

After convergence, samples that have reached true local maxima are marked as stuck: \[\mathrm{stuck}_i \gets \left( \|\nabla f(x_i)\|_\infty < \epsilon_{\mathrm{grad}} \right) \wedge \left( f(x_i) \geq f_{\max} + \log \delta \right)\] The gradient threshold \(\epsilon_{\mathrm{grad}}\) (default \(10^{-6}\)) ensures the sample has reached a stationary point. The likelihood threshold \(\log \delta\) (default \(\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.3 Phase 3: Deduplication↩︎

Stuck samples are deduplicated using the \(L_\infty\) metric (see Appendix 7): \[\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)\] The sample with the higher function value survives; its inverse Hessian estimate is preserved for width estimation. The number of removed duplicates \(K_{\mathrm{lost}}\) is recorded.

3.2.4 Phase 4: Reseeding (Repulse Monkey / Golden Rooster)↩︎

If \(K_{\mathrm{lost}} > 0\) and this is not the last oscillation, the population is maintained at size \(N\) by reseeding. Two strategies are used depending on exhaustion state.

3.2.4.1 Repulse Monkey.

When \(\geq 5\) unfrozen 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.

3.2.4.2 Golden Rooster.

When \(< 5\) unfrozen 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.5 Phase 5: Hands Like Clouds↩︎

If this is not the last oscillation, unfrozen samples take \(n_{\mathrm{cloud}}\) gradient ascent steps using a stochastically smoothed gradient: \[\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) \label{eq:smoothed-grad}\tag{1}\] Smoothing at scale \(\sigma\) blurs features smaller than \(\sigma\) while preserving global basin structure. For multi-scale problems (Rastrigin, Ackley), this allows samples to ascend toward global basins through local ripples. The smoothing scale \(\sigma\) is auto-estimated from sample spread: \(\sigma = c \cdot \mathrm{std}(\{x_i : \neg\mathrm{stuck}_i\})\) with \(c = 0.1\) by default, or supplied by the caller. Stuck samples are excluded from this phase.

3.2.6 Phase 6: Anti-Convergence↩︎

Unfrozen samples receive \(n_{\mathrm{anti}}\) momentum-based gradient descent steps: \[v_{t+1} = \mu \cdot v_t - \alpha \cdot \nabla f(x_t), \qquad x_{t+1} = x_t + v_{t+1} \label{eq:anticonverge}\tag{2}\] with momentum \(\mu = 0.9\) and 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 = n_{\mathrm{osc}}\)), so that the final convergence phase yields clean output peaks without further perturbation.

3.3 Design 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.4 Complete Algorithm↩︎

3.5 GPU Execution Model↩︎

Every operation in Algorithm 1 vectorizes over the sample index \(i\). Let \(P\) denote the number of GPU parallel processors (NVIDIA RTX 3080: \(P = 8704\) CUDA cores). For a population of size \(N\):

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

3.6 Hyperparameters↩︎

ChiSao has several hyperparameters. Table 1 lists defaults and sensitivity range (see Appendix 9 for full sensitivity analysis).

4 Theoretical Analysis↩︎

4.1 Log-Concave Case↩︎

4.2 Multimodal Case: Coverage Analysis↩︎

For non-log-concave \(f\) with \(K^*\) significant modes \(\{\theta^*_k\}_{k=1}^{K^*}\), we analyze the probability that all modes are discovered.

Let \(\mathcal{B}_k = \{\theta : \theta^*_k = \mathop{\mathrm{arg\,max}}_{\theta'} f(\theta'), \mathrm{L\text{-}BFGS}(\theta) \to \theta^*_k\}\) denote the basin of attraction of mode \(k\) under L-BFGS, with volume \(V_k = \mathrm{vol}(\mathcal{B}_k \cap \Theta)\) and total basin volume \(V_{\mathrm{basin}} = \sum_k V_k\).

4.3 GPU Complexity↩︎

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

4.4 Limitations of the Analysis↩︎

The convergence guarantee (Proposition 1) requires log-concavity, which fails for all genuinely multimodal functions. The coverage analysis (Proposition 2) 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.

5 Experiments↩︎

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). ChiSao is implemented in Python using CuPy [22] for GPU execution, with NumPy [23] fallback. All hyperparameters use the defaults in Table 1; no per-problem tuning was performed. Every benchmark function is evaluated value-only: ChiSao computes its gradients by finite differences, so every wall-clock time reported below is the derivative-free worst case [24]. Analytic or automatic-differentiation gradients, when available, only reduce it.

5.1 Benchmark Functions↩︎

We evaluate on all 42 standard test functions from the Simon Fraser University optimization benchmark suite [25]; see [26] for the wider catalogue and [27] for the CEC niching competition suite, spanning qualitatively distinct landscape types. Functions are organized into four groups (Table 2): scalable multimodal (Group A, tested across \(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 size \(N = 200\).

5.2 Metrics↩︎

We report two metrics. Mode recovery rate: fraction of trials in which the global optimum is located within a function-specific \(L_\infty\) tolerance (Table 2). Mean wall-clock time: total runtime in seconds averaged over 10 trials.

5.3 Baselines↩︎

We compare ChiSao against three standard baselines. Differential Evolution (DE): SciPy implementation [28] with default strategy. Basin-Hopping (BH): SciPy implementation with step-size adaptation. CMA-ES: cma package with multi-start restarts. ChiSao is evaluated with two seeding strategies: random (uniform random initialization) and carry_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]). All methods use matched function evaluation budgets.

5.4 Mode Recovery Results↩︎

Tables 3–5 report mode recovery rates for \(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.

5.4.0.1 Scalable multimodal functions (Group A).

The dominant result is consistent 100% recovery by both ChiSao seeders on Rastrigin, Levy, Styblinski-Tang, and Michalewicz through \(d = 64\). On Rastrigin, DE achieves 90% at \(d = 2\) but collapses entirely at \(d \geq 8\); BH and CMA-ES follow the same pattern. On Levy and Styblinski-Tang, DE and CMA-ES remain competitive through \(d = 8\) before collapsing, while both ChiSao seeders maintain perfect recovery. Michalewicz—with \(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).

The Ackley result reveals a structural difference between the two seeding strategies. Random seeding achieves 80% at \(d = 4\) and collapses entirely at \(d \geq 8\); carry_tiger maintains 100% from \(d = 2\) through \(d = 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% through \(d = 8\) before collapsing at \(d \geq 16\).

On Griewank, carry_tiger is the only method achieving 100% at \(d = 2\) and \(d = 4\) (random seeding: 10% and 0% respectively; DE: 60% and 10%; CMA-ES: 50% and 0%). From \(d \geq 8\) onward, both ChiSao seeders fully recover, and DE partially recovers at \(d \geq 32\) (80%). The carry_tiger advantage at low-\(d\) Griewank 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 for ChiSao. Both seeders fail at \(d \geq 16\) (reporting no peaks, denoted \(\dagger\)); carry_tiger achieves only 10% at \(d = 8\) and 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 threshold \(\log 0.1\). DE achieves 80% at \(d = 8\) before also collapsing at \(d \geq 16\); no method succeeds beyond \(d = 8\). This failure is outside ChiSao’s design envelope, discussed in Section 6.

5.4.0.2 Scalable bowl and valley functions (Group B).

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

On Rosenbrock, both seeders fail at \(d = 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 at \(d \geq 4\). DE and BH succeed at \(d = 2\); CMA-ES collapses from \(d = 8\).

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

Sum of Different Powers shows a subtler pattern: ChiSao random drops to 60% at \(d = 64\) while carry_tiger remains at 100%, confirming that structured ray coverage provides a meaningful advantage at the largest dimensions.

5.4.0.3 Fixed-2D multimodal functions (Group C).

Table 5 reports results for all 13 Group C functions. Both ChiSao seeders achieve 100% on all but two. Easom (20% random, 60% CT) presents an isolated basin covering roughly \(10^{-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 high \(d\); DE, BH, and CMA-ES all achieve 100%. Drop Wave is the only function where ChiSao outperforms a pure baseline: both seeders achieve 100% while DE reaches only 60% and CMA-ES 40%.

5.5 Wall-Clock Scaling↩︎

Table 6 reports mean wall-clock times for six representative functions covering the qualitatively distinct outcomes in the benchmark: functions where only ChiSao succeeds at high \(d\) (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\%\) recovery at that dimension; exact recovery rates are given in Tables 3–4.

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 across \(d = 2\) to \(d = 64\); all baselines require increasing time for 0% recovery from \(d \geq 8\) onward. On Styblinski-Tang, ChiSao (random) completes in 2.8–5.2 s across the full range while DE and CMA-ES collapse at \(d \geq 32\).

On Levy at \(d = 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. At \(d = 64\), ChiSao (random) completes in 23.8 s and carry_tiger in 29.2 s; BH requires 805.4 s (\({\sim}34\times\) slower), CMA-ES 46.5 s (\(2\times\)), and DE 36.1 s (\(1.5\times\)). The \(34\times\) BH–ChiSao ratio at \(d = 64\) represents 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: at \(d = 64\), all four methods achieving 100% require 50.2 s (ChiSao random), 99.0 s (CT), 312.5 s (DE), and 1946.4 s (BH), a \(39\times\) gap between ChiSao and 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=2\)), 7.0 s (\(d=8\)), 4.6 s (\(d=32\)), 23.1 s (\(d=64\)). The decrease from \(d=2\) to \(d=32\) reflects early oscillation-cycle termination when carry_tiger seeding places samples close to the global basin, reducing the number of exploration phases required. The \(d=64\) increase reflects additional L-BFGS iterations needed to converge in the high-dimensional flat landscape.

5.6 Ablation Study↩︎

Three results dominate the ablation (Table 7), with the \(d=8\) findings extended by additional ablations at \(d \in \{16, 32, 64\}\) under carry_tiger seeding.

5.6.0.1 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.

5.6.0.2 Full ChiSao is the most robust configuration across dimensions.

Across all tested dimensions (\(d \in \{8, 16, 32, 64\}\)), full ChiSao with 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 at \(d \geq 32\) as 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.

5.6.0.3 Anti-convergence effects are dimension-dependent.

At \(d=8\) under 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: at \(d=16\), Levy recovery under the Ph. 6-off configuration drops to 20%; at \(d=32\), Levy fails entirely while Schwefel recovers only 60%; at \(d=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 at \(d=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. The \(d=8\) ablation 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.7 Noise Robustness↩︎

We test ChiSao’s mode detection reliability when \(f\) is corrupted by additive Gaussian noise: \(\tilde{f}(\theta) = f(\theta) + \varepsilon\), \(\varepsilon \sim \mathcal{N}(0, \sigma_{\mathrm{noise}}^2)\). We use a 2-mode Gaussian mixture in \(d = 6\) with peaks at \(\pm 2\hat{e}_1\), mode separation \(4\sigma\).

These results are available from the SunBURST benchmark suite [29]; we reproduce them here for completeness.

The key result is 100% mode detection even at \(\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 \(\sigma = 0.5\). This robustness follows from ChiSao’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.

6 Discussion↩︎

The noise result is the one that matters (Table 8). 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.

6.0.0.1 Relation to SunBURST.

ChiSao was developed as the mode-discovery engine for SunBURST [29], 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.

6.0.0.2 Limitations.

ChiSao assumes 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 \(\sigma\) is auto-estimated but may require manual tuning on strongly heterogeneous landscapes. The anti-convergence momentum step size \(\alpha\) is sensitive on functions with very small or very large gradient magnitudes; an adaptive line-search would improve robustness. The theoretical analysis (Section 4) covers only the log-concave and first-oscillation-pass cases; a full convergence theory for the oscillation dynamics remains open.

6.0.0.3 Failure modes.

The quality-gate threshold (\(\delta = 0.1\), i.e., modes within one log-decade of the global maximum) is the primary source of ChiSao’s failures on this benchmark. On Schwefel, the deceptive landscape drives convergence to secondary peaks whose likelihood falls below threshold; on Trid at high \(d\) and 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 where ChiSao fails on these functions, at the cost of reporting spurious peaks on harder multimodal problems. ChiSao underperforms additionally on functions with exponentially many equal-height modes: the deduplication and reseeding logic does not scale to \(O(2^d)\) distinct modes.

6.0.0.4 Future work.

Since all reported results already use finite-difference gradients, analytic or automatic-differentiation gradients are a direct speedup lever: they remove the \(2dN\) finite-difference evaluations per L-BFGS step, though the GPU already absorbs them in \(O(1)\) wall-clock while \(2dN < 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 and the empirical multi-modal results.

Data Availability↩︎

All benchmark scripts, result files, and analysis code are available at https://github.com/beastraban/chisao.

Code Availability↩︎

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

Competing Interests↩︎

The author declares no competing interests.

7 The \(L_\infty\) Metric for Deduplication↩︎

The \(L_\infty\) (Chebyshev) distance \(\|x - y\|_\infty = \max_i |x_i - y_i|\) has two advantages over \(L_2\) for high-dimensional deduplication. First, it is dimension-independent in expectation: for uniform random points in \([0,1]^d\), \(\mathbb{E}[\|x-y\|_\infty] \approx d/(d+1) \to 1\), while \(L_2\) distance diverges as \(\sqrt{d/6}\). A fixed threshold \(\epsilon_{\mathrm{dup}}\) in \(L_\infty\) therefore has consistent semantics across dimensions. Second, the \(L_\infty\) ball \(\{y : \|x-y\|_\infty \leq r\}\) is a hypercube of side \(2r\), enabling efficient spatial hashing for \(O(N)\) deduplication rather than \(O(N^2)\) pairwise comparison.

8 Hands Like Clouds: Smoothed Gradient Analysis↩︎

The smoothed function \(f_\sigma(\theta) = \mathbb{E}_{z \sim \mathcal{N}(0,I_d)}[f(\theta + \sigma z)]\) satisfies: \[\nabla f_\sigma(\theta) = \mathbb{E}_{z \sim \mathcal{N}(0,I_d)}[\nabla f(\theta + \sigma z)]\] by Leibniz’s rule (assuming \(f\) is differentiable a.e.). The Monte Carlo estimator 1 is unbiased. For functions with multi-scale structure \(f = f_{\mathrm{global}} + f_{\mathrm{local}}\) where \(f_{\mathrm{local}}\) has characteristic scale \(\lambda < \sigma\), the smoothed gradient \(\nabla f_\sigma\) eliminates \(f_{\mathrm{local}}\) contributions (their integral over \(\mathcal{N}(0,\sigma^2 I)\) vanishes by symmetry for zero-mean local oscillations), revealing the global basin structure.

[18] showed that for convex \(f\), gradient descent on \(f_\sigma\) converges to within \(O(\sigma^2 L)\) of the optimum of \(f\), where \(L\) is the Lipschitz constant of \(\nabla f\). For non-convex \(f\), the smoothing provides no convergence guarantee but empirically helps samples navigate past local barriers whose scale is smaller than \(\sigma\).

9 Hyperparameter Sensitivity↩︎

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

Three observations follow. (i) On the six landscapes that ChiSao solves at 100% with the defaults, recovery is invariant to every perturbation in the sweep with one exception. (ii) That exception is \(n_{\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 default \(n_{\mathrm{anti}} = 5\) is 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: Schwefel’s failure is a property of the quality-gate threshold \(\delta = 0.1\) interacting with the deceptive landscape, not of the search-dynamics parameters that the sweep covers. A sensitivity sweep over \(\delta\) would close this loop, but \(\delta\) is not currently exposed as a hyperparameter in the sticky_hands interface.

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

10 Benchmark Function Definitions↩︎

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

10.0.0.1 Rastrigin.

\(f(\theta) = -\left[Ad + \sum_{i=1}^d (\theta_i^2 - A\cos(2\pi\theta_i))\right]\), \(A=10\), \(\theta \in [-5.12, 5.12]^d\). Global maximum at \(\theta = 0\), \(f(0) = 0\).

10.0.0.2 Ackley.

\(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=20\), \(b=0.2\), \(c=2\pi\), \(\theta \in [-32.768, 32.768]^d\).

10.0.0.3 Schwefel.

\(f(\theta) = -\left[418.9829d - \sum_{i=1}^d \theta_i \sin(\sqrt{|\theta_i|})\right]\), \(\theta \in [-500, 500]^d\). Global maximum near \(\theta_i = 420.97\).

10.0.0.4 Styblinski-Tang.

\(f(\theta) = -\frac{1}{2}\sum_{i=1}^d (\theta_i^4 - 16\theta_i^2 + 5\theta_i)\), \(\theta \in [-5, 5]^d\).

10.0.0.5 Michalewicz.

\(f(\theta) = \sum_{i=1}^d \sin(\theta_i)\sin^{2m}(i\theta_i^2/\pi)\), \(m=10\), \(\theta \in [0, \pi]^d\).

10.0.0.6 Trid.

\(f(\theta) = -\left[\sum_{i=1}^d (\theta_i - 1)^2 - \sum_{i=2}^d \theta_i \theta_{i-1}\right]\), \(\theta \in [-d^2, d^2]^d\). Global maximum at \(\theta_i^* = i(d+1-i)\) with \(f^* = d(d+4)(d-1)/6\). Unimodal but increasingly ill-conditioned with \(d\): Hessian eigenvalues span many orders of magnitude and the optimum coordinates \(\theta_i^*\) range from \(d\) to \(d(d+1)^2/4\).

10.0.0.7 Bukin N.6.

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

10.0.0.8 Gaussian mixture.

\(f(\theta) = \log \sum_{k=1}^K w_k \mathcal{N}(\theta; \mu_k, \sigma^2 I)\) with \(w_k = 1/K\), \(\mu_k\) random in \(\Theta\), \(\sigma = 1\), \(\Theta = [-10, 10]^d\).

References↩︎

---

1. Wing 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.↩︎