Combinatorial and continuous optimization problems with NP-hard characteristics remain among the most challenging computational problems across scientific and engineering fields . Classical optimization techniques often struggle with the fundamental trade-off between exploration and exploitation, where global search methods risk computational inefficiency while local search heuristics frequently converge prematurely to suboptimal solutions . Markov Chain Monte Carlo (MCMC) methods provide powerful mechanisms for global exploration through probabilistic sampling of the search space , particularly effective in high-dimensional and multimodal landscapes. However, MCMC approaches alone lack the aggressive local refinement needed for rapid convergence to higher-quality solutions. Conversely, greedy local search methods excel at exploitation of promising regions but offer no systematic mechanism for escaping local optima. Simulated Annealing (SA) addresses this problem through temperature-controlled stochastic acceptance, though it requires careful parameter tuning and may fail to escape deep local minima without adaptive reheating techniques .

Existing state-of-the-art optimizers each address different aspects of this problem. Bayesian Optimization excels in low-dimensional smooth landscapes but scales poorly. Covariance Matrix Adaptation Evolution Strategy (CMA-ES) provides robust derivative-free optimization but incurs significant computational overhead. Accelerated Particle Swarm Optimization (APSO) offers good exploration-exploitation balance but lacks mechanisms to avoid revisiting poor regions. Genetic Algorithms introduce population-based evolutionary computation with selection, crossover, and mutation operators, while Tabu Search employs memory structures to prevent cycling and encourage exploration. No single classical approach effectively combines global exploration, local exploitation, and adaptive escape mechanisms while maintaining computational efficiency across diverse problem classes.

We present Yukthi Opus (YO), a three-layer hybrid metaheuristic optimizer that systematically integrates MCMC-based global exploration, greedy local search, and adaptive simulated annealing with reheating, following the memetic algorithm paradigm of combining population-based and local search methods. YO addresses several critical gaps: preventing premature convergence through structured burn-in exploration, avoiding computational waste via blacklist mechanisms that prevent revisiting poor regions, escaping local minima through adaptive temperature reheating, and maintaining solution robustness through multi-chain parallel execution with post-burnin selection.

Our key contributions include introducing a novel three-layer hybrid design that combines MCMC , greedy search, and SA with adaptive reheating in a principled structure allowing explicit control over evaluation budgets, implementing a spatial blacklist system that prevents repeated evaluation of demonstrably poor parameter regions, demonstrating through experiments that parallel chain execution with post-burnin selection improves solution stability and reduces variance compared to single-chain approaches, conducting comprehensive ablation studies on the Rastrigin 5D function to quantify the individual contributions of MCMC, greedy search, SA, blacklisting, and multichain execution, evaluating YO on the Traveling Salesman Problem with 50-200 cities across multiple random seeds, and benchmarking YO against state-of-the-art methods including CMA-ES , Bayesian Optimization , and APSO on the challenging Rosenbrock 5D function.

Overview

The Yukthi Opus (YO) Hybrid Optimizer is a three-layer metaheuristic designed for NP-hard optimization by integrating complementary search strategies . It explicitly controls the evaluation budget, allocating resources between an exploratory MCMC burn-in phase and a hybrid exploitation phase, making it well suited for expensive black-box objectives. The multi-chain structure improves robustness and reduces sensitivity to initialization.

Initialization defines the search bounds, total budget, parallel chains, random starting points, blacklist, and simulated annealing (SA) temperature schedule. Phase 1 performs MCMC burn-in for global exploration, independently sampling and accepting candidates per chain using the Metropolis criterion while tracking the best solutions and optionally blacklisting poor regions. The top-performing samples seed Phase 2, which combines MCMC proposals, greedy local refinement, and SA-based acceptance with cooling and optional reheating under stagnation. This design balances global exploration with aggressive local convergence.

Finally, post-processing aggregates results across chains, selects the best global solution, computes performance metrics, and returns the optimal solution along with the optimization trace.

Algorithm Architecture and Flow

Operational Workflow

Figure 1 illustrates YO’s operational workflow, showing the two-phase architecture with adaptive mechanisms.

Algorithmic Pseudocode

Core Components and Design Rationale

YO integrates six complementary components, each addressing specific optimization challenges. Table 1 summarizes their descriptions and contributions.

Classification as Hybrid Metaheuristic

YO is classified as a hybrid metaheuristic because it systematically combines multiple distinct optimization methods: stochastic global search (MCMC), deterministic local search (greedy), and adaptive stochastic acceptance (SA with reheating). Unlike single-strategy metaheuristics that rely on one method such as pure genetic algorithms or pure simulated annealing , hybrid metaheuristics leverage complementary strengths of different approaches to address the exploration-exploitation trade-off more effectively. The three-layer design ensures that exploration and exploitation occur in structured phases with explicit resource allocation, rather than competing for evaluations in an uncoordinated manner. The blacklist and post-burnin selection further enhance efficiency by directing computational budget toward promising regions. This principled integration of complementary components, explicit budget control, and adaptive escape mechanisms distinguish YO from both classical single-strategy metaheuristics and simple ensemble approaches.

Experimental Results

We evaluate YO on three challenging NP-hard benchmarks: the Rastrigin 5D function with comprehensive ablation studies, the Traveling Salesman Problem with 50-200 cities, and the Rosenbrock 5D function with state-of-the-art comparisons. All results presented are taken directly from experimental data without modification.

Rastrigin 5D Function: Ablation Studies

The Rastrigin function is ideal for ablation studies due to its highly multimodal landscape with numerous regularly distributed local minima separated by the same barrier height. The 5D version provides sufficient complexity to test whether YO’s modules actually contribute to improved convergence while remaining computationally tractable for repeated runs. We systematically disable individual components to isolate their contributions to solution quality, convergence speed, and stability.

Experimental Setup

The problem uses dimensionality $`D=5`$, search space bounds $`[-5.12, 5.12]^5`$, evaluation budget of 150 evaluations per run, and 30 runs per variant for statistical significance. The test function is an expensive multi-modal function combining Rastrigin (multiple local minima), Rosenbrock (narrow valley), Sphere (convex bowl), and sin plus exponential terms, with delay of 0.01 seconds per evaluation to simulate expensive black-box functions. We test six ablation variants: A0_Full_YO as complete baseline, A1_No_MCMC removing MCMC exploration phase, A2_No_Greedy removing greedy local search refinement, A3_No_SA removing simulated annealing acceptance control, A4_No_Blacklist disabling blacklist mechanism, and A5_Single_Chain using only one chain instead of multiple parallel chains. All other parameters are held constant to isolate effects.

Quantitative Results

Table 2 presents ablation study results with statistical significance tested using two-sample t-tests.

Note: CV = Coefficient of Variation (Std/Mean). Statistical significance: *** = $`p < 0.01`$, * = $`p < 0.05`$.

Figure 2 shows distribution of results across all variants.

Component-wise Analysis

Removing MCMC exploration causes the most severe degradation with 36% worse solution quality and 26% reduction in stability (CV increases from 0.331 to 0.417). This demonstrates MCMC is critical for global exploration in multimodal landscapes. Without MCMC, the optimizer relies solely on greedy search from random initialization, frequently converging to poor local minima. While removing MCMC provides 32% runtime speedup, this is massively outweighed by quality loss.

Removing greedy local search causes 30% degradation in mean solution quality. Interestingly, stability improves (CV: 0.207 vs 0.331), suggesting greedy search introduces variance by aggressively exploiting diverse local basins discovered by MCMC. Without greedy refinement, the optimizer relies on stochastic MCMC proposals and SA acceptance, converging more slowly but more uniformly across runs. This demonstrates greedy’s critical role in rapid local refinement.

Removing Simulated Annealing causes 25% degradation and significantly reduced stability (CV: 0.438 vs 0.331). Without SA’s temperature-controlled acceptance mechanism, the optimizer cannot effectively balance exploration and exploitation or escape from suboptimal basins. The high CV indicates highly inconsistent performance as the algorithm occasionally gets trapped in deep local minima without probabilistic uphill moves.

Removing the blacklist mechanism has no measurable impact (identical mean, std, CV to baseline). This suggests the Rastrigin 5D landscape, while highly multimodal, does not contain large spatially contiguous poor regions that would be repeatedly revisited. The blacklist provides value in problems with pathological regions but not in uniformly distributed multimodal landscapes.

Using a single chain paradoxically improves mean performance by 30% but drastically reduces stability (CV increases from 0.331 to 0.734, more than doubling). This high variance indicates single-chain occasionally finds excellent solutions through fortunate initialization but frequently performs poorly when initialized in bad regions. The multi-chain design trades slight degradation in best-case performance for substantial variance reduction, a worthwhile trade-off for production applications requiring consistent performance.

The ablation study reveals that MCMC and Greedy are critical components forming the core of YO’s effectiveness (removing either causes statistically significant degradation of 30-36%). SA and multi-chain execution primarily improve robustness and consistency rather than best-case performance (removal increases variance by 32-122%). The blacklist provides context-dependent value. The full YO pipeline provides the best balance with no obviously redundant components.

Traveling Salesman Problem (TSP)

The TSP benchmark evaluates YO on Euclidean instances seeking the shortest tour visiting all cities exactly once and returning to start, a canonical NP-hard combinatorial optimization problem. We test three problem sizes (50, 100, 200 cities) with evaluation budgets scaled accordingly (20,000 for $`N=50`$, 50,000 for $`N=100`$, 100,000 for $`N=200`$). Three random seeds (42, 101, 202) are used for each problem size to assess stability.

Baseline Methods

YO is compared against four established TSP heuristics: Simulated Annealing (classical SA with temperature annealing), Genetic Algorithm (population-based evolutionary approach), 2-opt Restart (deterministic local search with random restarts), and Random Search (uniform random sampling baseline). YO uses MCMC for exploration, greedy 2-opt for exploitation, blacklist to avoid poor tours, adaptive reheating to escape local minima, and multi-chain approach for robustness.

Solution Visualization

Table [tab:tsp_aggregate] presents aggregate statistics across three seeds for each problem size.

Analysis

TSP results demonstrate problem-size-dependent performance. For small instances ($`N=50`$), YO achieves tour quality statistically equivalent to 2-opt Restart (5404.10 vs 5403.49 mean, difference $`<0.02\%`$) but requires significantly higher runtime (25.42s vs 11.67s, 2.2$`\times`$ slower). The overhead from MCMC exploration and multi-chain execution provides negligible benefit for small TSP instances where deterministic local search with restarts suffices.

For medium instances ($`N=100`$), YO shows clear advantages finding tours 0.5% shorter than 2-opt (7620.21 vs 7657.17 mean) with comparable standard deviation (331.02 vs 349.31). Runtime increases to 193.91s compared to 2-opt’s 124.82s (1.55$`\times`$ slower), but YO substantially outperforms SA (9.2% improvement) and GA (61.0% improvement).

For large instances ($`N=200`$), YO achieves best-in-class performance finding tours 1.8% shorter than 2-opt (10715.36 vs 10914.85 mean) and 11.3% shorter than SA. Notably, YO exhibits lower variance (184.97) than 2-opt (194.51), indicating improved robustness as problem complexity increases. Runtime remains substantial at 1577.99s but the 60% runtime overhead relative to 2-opt (1.60$`\times`$ slower) is justified by consistent solution quality improvements.

The blacklist mechanism was triggered minimally (1-2 regions per run), suggesting TSP landscapes do not contain large spatially clustered poor regions under Euclidean distance metrics. YO’s strength on larger TSP instances stems from effective exploration of the exponentially growing solution space ($`|S| = (n-1)!/2`$), where simple restarts of deterministic local search become less likely to discover high-quality basins of attraction. The multi-chain design provides robust coverage of diverse tour topologies, while MCMC proposals enable transitions between local optima that pure 2-opt cannot traverse.

Rosenbrock 5D Function: State-of-the-Art Comparison

The Rosenbrock function, also known as the banana-shaped valley function, provides a challenging test case for gradient-free optimizers due to its narrow curved valley structure. We compare YO against three state-of-the-art optimizers: CMA-ES (Covariance Matrix Adaptation Evolution Strategy), BayesOpt (Bayesian Optimization with Gaussian Process surrogate), and APSO (Accelerated Particle Swarm Optimization). The narrow curved valley challenges exploration methods while favoring gradient-aware approaches, providing a demanding test case.

Experimental Setup

The problem uses dimensionality $`D=5`$, search space $`[-5.0, 10.0]^5`$ with global minimum $`f(1,1,1,1,1) = 0`$, evaluation budget of 150 evaluations per run, and 30 runs per optimizer for statistical analysis executed in parallel as independent runs.

Results Visualization

Figure 9 shows distribution of final best values across 30 runs for all optimizers.

Quantitative Performance

Table 6 presents quantitative results.

YO Hybrid achieves mean runtime 0.061s $`\pm`$ 0.012s (fastest among all optimizers), best solution 24.8463, median 602.2053, and coefficient of variation 1.458 (high variability). Table 7 presents statistical significance tests.

Performance Analysis

BayesOpt achieves by far the best solution quality (176.46 mean), significantly outperforming YO (1297.91, 7.4$`\times`$ worse), APSO (1489.69), and CMA-ES (50516.96). However, YO ranks second in solution quality while being fastest in runtime (0.061s mean), nearly 2$`\times`$ faster than competitors. YO demonstrates fastest runtime, excellent speed-accuracy trade-off for rapid optimization under tight time constraints, competitive best-case performance (minimum 24.85 close to BayesOpt 19.73), and robustness across initializations through multi-chain exploration.

However, YO exhibits weaknesses including high variability (standard deviation exceeding mean, CV=1.458) indicating inconsistent convergence, inability to match gradient-aware methods as BayesOpt’s Gaussian Process surrogate effectively exploits smooth valley structure, challenges with narrow valleys where random proposals struggle to follow curved valley while BayesOpt’s acquisition function naturally follows gradients, and greedy search limitations without gradient information making local refinement inefficient in narrow valleys.

Statistical tests confirm BayesOpt’s dominance is highly significant ($`p < 0.003`$) with large effect size (Cohen’s d = 0.837), indicating YO should not be first choice for Rosenbrock-like problems where practitioners should use BayesOpt instead. However, YO’s 2$`\times`$ speed advantage may be valuable when approximate solutions are acceptable. YO shines in scenarios with limited evaluation budgets (150-1000) where speed is critical, problems where objective evaluation is expensive relative to optimizer overhead, cases needing diverse exploration over absolute precision, and black-box optimization without available gradient or smoothness information. YO struggles with smooth low-dimensional problems where BayesOpt can build accurate GP surrogates, narrow valley structures where gradient-aware methods navigate efficiently, problems with strong local structure where covariance information or gradients help, and when consistency is critical as YO’s high CV (1.458) versus BayesOpt’s moderate CV (0.69) is problematic.

The overall assessment is that YO provides excellent runtime efficiency but cannot match BayesOpt’s solution quality on smooth, low-dimensional problems like Rosenbrock 5D. The narrow curved valley structure strongly favors gradient-aware methods like BayesOpt which can build surrogate models capturing the valley’s geometry. YO’s value proposition is the speed-accuracy trade-off for rapid optimization under constrained evaluation budgets, not absolute solution quality on smooth problems.

Discussion

Strengths of YO

YO demonstrates effective exploration-exploitation balance through its three-layer design systematically combining MCMC global exploration with greedy local exploitation, mediated by SA acceptance control. Ablation studies on Rastrigin 5D prove both components are critical for solution quality on multimodal problems, as removing either causes 30-36% performance degradation. The blacklist mechanism provides evaluation efficiency by preventing wasted evaluations on known poor regions, providing value on TSP by avoiding pathological subcycles. Multi-chain robustness substantially reduces solution variance, demonstrated quantitatively by Rastrigin ablation study showing 55% variance reduction (CV: 0.331 for multi-chain vs 0.734 for single-chain), valuable for production applications requiring consistent performance across runs. Adaptive escape mechanisms through reheating enable structured escape from local minima without manual intervention, with TSP and Rastrigin results showing clear convergence improvements following reheating events. YO exhibits domain generality, performing competently across diverse problem classes (combinatorial TSP, continuous multimodal functions) without problem-specific modifications beyond local search operators, making it reliable when problem structure is unknown. Explicit budget control distinguishes YO from population-based methods with unclear stopping criteria, providing explicit evaluation budget allocation between exploration (burn-in) and exploitation (hybrid optimization) phases, making it predictable for resource-constrained scenarios.

Weaknesses of YO

Computational overhead is substantial compared to simpler methods, with TSP $`N=50`$ showing 2.2$`\times`$ slowdown versus 2-opt, justified only when solution quality improvements compensate for increased runtime. Small problem inefficiency is evident, as for small TSP instances ($`N=50`$, only 0.02% improvement over 2-opt) and smooth landscapes (Rosenbrock), YO’s complexity provides negligible benefit over simple heuristics with overhead dominating on problems easily solvable by deterministic local search. Gradient-free limitations prevent YO from exploiting gradient information available in smooth problems, as Rosenbrock 5D results show BayesOpt significantly outperforms YO (176.46 vs 1297.91 mean, 7.4$`\times`$ better) by building gradient-aware surrogate models, making YO best suited for black-box optimization where gradients are unavailable or unreliable. High variance on some problems is demonstrated by Rosenbrock results showing very high coefficient of variation (1.458, with std exceeding mean) indicating inconsistent convergence, and while multi-chain execution mitigates this issue compared to single-chain (Rastrigin ablation: CV 0.331 vs 0.734), it does not eliminate the problem on narrow-valley landscapes.

Performance Characterization

Table 8 summarizes when YO shines versus when simpler alternatives suffice.

Key Design Insights

The Rastrigin ablation study provides quantitative evidence that multi-chain execution is critical for robustness through reducing coefficient of variation by 55% (from 0.734 to 0.331), providing consistent performance across runs, reducing sensitivity to initialization, and increasing probability of finding good basins through diverse exploration. This robustness stems from diverse initialization allowing different chains to explore different regions reducing dependence on lucky starting points, multiple independent searches increasing probability that at least one chain finds high-quality basin, best-of-N selection providing implicit variance reduction through selecting best result across chains, and portfolio effect ensuring poor performance of some chains compensated by good performance of others. For production applications requiring consistent performance (automated parameter tuning, deployment in critical systems), multi-chain overhead is justified by variance reduction.

Adaptive reheating addresses a fundamental limitation of classical simulated annealing where once temperature drops too low, the algorithm cannot escape deep local minima without external intervention. YO’s reheating mechanism monitors stagnation automatically with no manual tuning needed, temporarily increases acceptance probability to enable probabilistic uphill moves to escape local traps, enables structured recovery from premature convergence without restarting optimization, and provides empirically visible benefit as shown in TSP and Rastrigin results with clear convergence jumps following reheating events. The reheating mechanism is particularly valuable in problems where local minima have varying depths (some easy to escape, some difficult), search may stagnate mid-optimization despite remaining evaluation budget, and temperature cooling rate is difficult to tune a priori. Without reheating, classical SA requires either very slow cooling schedules (wasting evaluations) or risks premature convergence, with adaptive reheating providing middle ground.

Conclusion

We have presented Yukthi Opus (YO), a three-layer hybrid metaheuristic optimizer that systematically combines MCMC exploration, greedy local search, and adaptive simulated annealing with reheating. Comprehensive benchmarking across three diverse NP-hard problems (Rastrigin 5D with ablation studies, TSP with 50-200 cities, Rosenbrock 5D with state-of-the-art comparisons) establishes YO’s performance characteristics across problem classes.

YO accomplishes competitive performance on large complex problems (TSP $`N \geq 100`$ shows 0.5-1.8% improvement over 2-opt with lower variance), explicit control over evaluation budgets through principled exploration-exploitation allocation via structured burn-in and hybrid optimization phases, enhanced robustness through multi-chain execution (55% variance reduction, CV: 0.331 vs 0.734 compared to single-chain approach), structured escape from local minima through adaptive reheating mechanism enabling recovery from stagnation as demonstrated by visible convergence improvements, and domain generality being effective across combinatorial, continuous problems without domain-specific modifications beyond local search operators.

YO ranks second on Rosenbrock 5D but with fastest runtime offering excellent speed-accuracy trade-off, outperforms classical heuristics (SA, GA) on large problems but incurs overhead on small instances, and is complementary to specialized optimizers where practitioners should use BayesOpt for smooth low-dimensional problems and use YO for multimodal black-box optimization. YO is best suited for 100-1000 evaluation budgets, expensive objective functions where evaluation count is primary constraint, black-box problems without gradient information, multimodal landscapes with numerous local minima, and applications requiring consistent performance across runs. YO is not recommended for smooth low-dimensional problems (use BayesOpt), small instances solvable by simple heuristics, problems with available gradient information, and extreme runtime constraints where simplicity preferred.

The ablation studies demonstrate that MCMC and greedy components are critical for solution quality (removing either causes 30-36% degradation), while SA and multi-chain execution primarily improve robustness (reducing CV by 32-55%). The blacklist method provides value on problems with spatially clustered poor regions but has minimal impact on uniformly multimodal landscapes. Future work should focus on reducing computational overhead through surrogate models and multi-fidelity approaches, extending to constrained and multi-objective problems, developing theoretical convergence guarantees, and creating problem-specific templates and automated hyperparameter tuning. The empirical results establish YO as a viable optimizer for a specific niche: expensive black-box optimization of complex multimodal problems where robustness and evaluation efficiency are critical.

Computational Environment and Reproducibility

All experiments in this study were executed on Google Colaboratory (Colab) to ensure accessible, uniform, and reproducible evaluation across all benchmarks. Runs were performed on standard Colab runtimes equipped with 12 GB RAM and, when available, NVIDIA T4/P100-class GPUs. These resources were sufficient for all experiments, including large-scale TSP evaluations, N-body simulations, and multimodal function benchmarks. Python 3.8+ and commonly used scientific computing libraries (NumPy, SciPy, matplotlib, scikit-learn, and problem-specific packages) were used throughout.

To support full reproducibility, the complete set of benchmark implementations, optimizer code, experiment configurations, and plotting utilities is publicly available in the project repository:

The repository includes the exact scripts and notebooks used to generate all figures, tables, and results in this manuscript. Each experiment folder contains self-contained code that automatically installs required packages and executes the corresponding benchmark pipeline, ensuring that reproduction requires no additional setup beyond running the provided scripts in a Colab environment.

Acknowledgments

The authors thank the Department of ECE and Department of Science and Humanities at PES University for their support of this research.