arXiv CSopen access

BEAM: Bi-level Memory-adaptive Algorithmic Evolution for LLM-Powered Heuristic Design

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

1

BEAM: Bi-level Memory-adaptive Algorithmic Evolution for LLM-Powered Heuristic Design

arXiv:2604.12898v1 [cs.AI] 14 Apr 2026

Chuyang Xiang, Yichen Wei, Jiale Ma, Handing Wang, Senior Member, IEEE, Junchi Yan, Senior Member, IEEE

Abstract—Large Language Model-based Hyper Heuristic (LHH) has recently emerged as an efficient way for automatic heuristic design. However, most existing LHHs just perform well in optimizing a single function within a pre-defined solver. Their single-layer evolution makes them not effective enough to write a competent complete solver. While some variants incorporate hyperparameter tuning or attempt to generate complex code through iterative local modifications, they still lack a high-level algorithmic modeling, leading to limited exploration efficiency. To address this, we reformulate heuristic design as a Bi-level Optimization problem and propose BEAM (Bi-level Memory-adaptive Algorithmic Evolution). BEAM’s exterior layer evolves highlevel algorithmic structures with function placeholders through genetic algorithm (GA), while the interior layer realizes these placeholders via Monte Carlo Tree Search (MCTS). We further introduce an Adaptive Memory module to facilitate complex code generation. To support the evaluation for complex code generation, we point out the limitations of starting LHHs from scratch or from code templates and introduce a Knowledge Augmentation (KA) Pipeline. Experimental results on several optimization problems demonstrate that BEAM significantly outperforms existing LHHs, notably reducing the optimality gap by 37.84% on aggregate in CVRP hybrid algorithm design. BEAM also designs a heuristic that outperforms SOTA Maximum Independent Set (MIS) solver KaMIS. Index Terms—Large Language Model, Heuristic Design, Metaheuristic, Monte Carlo Tree Search, Knowledge Augmentation.

I. I NTRODUCTION Heuristics are crucial for solving complex optimization problems, but manual design is laborious and biased [1]. Automatic Heuristic Design (AHD) emerged to mitigate this issue, with Hyper-Heuristics (HH) [2] automating parameter tuning and components combination—though inflexible. The rise of Large Language Model (LLM)-based code generation [3], [4] has opened up a new gate for AHD, yet general prompting strategies [5], [6] and general LLM agents fall short for this feedback-intensive task. Language Hyper-Heuristics (LHH) advances AHD by integrating LLM into frameworks such as the Genetic Algorithm (GA) [7]. In this line of work, heuristics are treated as individuals, and LLMs are used to improve these individuals iteratively guided by specialized prompts. However, existing LHHs only perform well in generating a single function of an algorithm [8] instead of entire ones, still demanding manual framework design. This reflects two fundamental limitations of these approaches: 1) Structural and Prompting Strategy Deficiencies: Most LHHs are single-layered, treating algorithms as single individuals. Chuyang Xiang and Yichen Wei contributed equally to this work.

Fig. 1: LHH Usage: 1) Single function design (used in [7], [11], [12]) for a given algorithm structure. 2) Hybrid algorithm design (proposed): designing a whole solver with given heuristic components. 3) Entire algorithm design: designing algorithms from scratch.

When complex requirements are given, their output codes remain simplistic, and the heuristics often fail to evolve after a few generations. While some variants attempt to generate complex code through iterative local modifications [9], they still lack a high-level algorithmic modeling. These frameworks may also degrade into random search as LLM cannot discern performance causality [7] when faced with complex codes. 2) Absent or Deficient Knowledge Augmentation: Existing approaches either let LLM design algorithms entirely from scratch—providing little or even zero textual external knowledge [7] or warm-start them with template functions [10] which demands significant manual intervention to design sufficiently diverse templates. To overcome these intertwined challenges, we argue that the automated generation of complex heuristics must align more closely with human algorithmic design principles. Human experts rarely construct sophisticated solvers as monolithic entities from scratch; instead, they decompose the problem into high-level structural planning (i.e., the algorithmic framework) and low-level component realization, frequently reusing and recombining established strategies [13]. Consequently, we advocate for a paradigm shift in LHHs that reformulates algorithm design as a bi-level optimization problem, allowing specialized search strategies to independently conquer framework evolution and function implementation. Furthermore, to prevent the LLM from conducting blind code exploration, this bi-level search must be firmly grounded in structured external knowledge and a repository of reusable heuristic components. This approach effectively bridges the generative flexibility of modern LLMs with the robust algorithmic recombination principles of traditional HH [14].

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

2

Fig. 2: Pipeline of our BEAM. The exterior layer designs the heuristic (code) structure (the heuristic()) via genetic evolution; the interior layer designs the functions (the func i()s) by MCTS. Evaluation happens every time when designing a new function to ensure the function quality. After the entire genetic evolution, the best heuristic (code structure + functions) are printed out, and the functions (the func i()s) are stored into the adaptive memory for future code structures to call directly. As shown in our experiments (Table V), BEAM has better performance over Google’s AlphaEvolve [9].

To address these limitations and realize this philosophy, we make the following contributions: • We propose BEAM (Bi-level Memory-adaptive Algorithmic Evolution) as shown in Section III, which reformulates AHD as a Bi-level Optimization problem [15], decomposing it into high-level structure generation (via GA) and low-level function realization (via MCTS). It’s further enhanced by an Adaptive Memory mechanism, enabling LLM to directly call elite low-level functions from previous generations. • We introduce a general Knowledge Augmentation (KA) pipeline (in Section IV-B) where LLM builds 2 datasets: a HeuBase of callable functions and a textbased KnoBase after retrieving external knowledge. We also construct part of HeuBase to incorporate cuttingedge heuristic components unavailable via pip, bridging traditional HH idea [14] with modern LLM capabilities. • We integrated the KA pipeline to our BEAM and baseline LHHs and tested them on desigining complete solvers for a series of combinatorial and continuous optimization problems. BEAM demonstrates significant performance improvements over existing LHHs and even surpassing SOTA solvers in MIS. In CVRP hybrid algorithm design, it delivers 37.84% aggregate advancement across all benchmarks. II. R ELATED W ORK Prompt Engineering for LLM Coding With the rapid progress of LLM in code generation [3], prompt engineering has emerged as a simple yet effective enhancement approach [4], [5]. Methods like CoT [6] and ToT [16] help

structure reasoning, while modular-inspired prompting (e.g. sketch-refine) improves control flow planning [17]. However, these general-purpose strategies often lack real-time feedback, limiting their effectiveness for black-box optimization tasks where high-quality heuristic generation is crucial [18]. LLM For Optimization Problems (LLM4OP). While LLM are limited in directly solving complex optimization problems [19], they excel at problem modeling and code generation. LLM4OP generally falls into two main categories: 1) Solver Assistance: LLM translate natural language into formal problem formulations and work with Neural CO solvers [20], directly generate solutions [21], [22] or solver-ready code [23], [24] using techniques like RAG [25], [26] or work with LSTM [27] to choose algorithms [28]. 2) Automatic Heuristic Design (AHD): LLM aid in designing new algorithms or heuristic components. Our proposed framework belongs to the AHD category. Language Hyper Heuristics (LHH). Early LHH attempts built a program database and let LLM iteratively refine them [29]. Later on, researchers were inspired by Hyper Heuristics (HH) and employed LLM as genetic operators [30] to evolve new heuristics. Evolution of Heuristics (EoH) used five fixed prompts to do this [11], [31], focusing on designing priority functions within predefined frameworks such as Guided Local Search (GLS), a metaheuristic that uses penalty-based guidance to escape local optima. Reflective Evolution (ReEvo) was an improved structure with redesigned prompts and reflection mechanism [7]. LLaMEA introduced a more statistically-sound evaluation method [8]. However, population-based search methods [32] often struggle to fully exploit the strengths of individual heuristics [33]. Recent work

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

3

TABLE I: Language Hyper Heuristic (LHH) comparison. *: Simple evolution means simply updating all heuristics using different prompts, †: AlphaEvolve outputs formatted instructions on revising code templates. Methods

Individual Strategy

Search Method

Calibration

FunSearch [29] EoH [11] EvoCAF [40] HSEvo [41] LLaMEA-HPO [42] MCTS-AHD [33] ReEvo [7] PoH [12] CPro1 [43] AlphaEvolve [9] BEAM (Ours)

One-shot One-shot One-shot One-shot One-shot One-shot Two-step (text reflection & codes) Two-step (text plans & codes) Two-step (text outline & codes) One-shot† Bi-level (algorithm structure & function realization)

Random sampling Simple evolution∗ Simple evolution GA GA MCTS GA MCTS Random sampling Random sampling CMA-ES

/ LLM LLM HS SMAC3 / / / Optuna /

Methods

KA Type

FunSearch [29] EoH [11] EvoCAF [40] HSEvo [41] LLaMEA-HPO [42] MCTS-AHD [33]

Templates Templates / Text / /

ReEvo [7]

PoH [12] CPro1 [43] AlphaEvolve [9] BEAM (Ours)

Benchmark Problems

Type

BPP BPP, TSP(GLS), FSSP(GLS) CAF BPO, TSP(GLS), OP BBOB TSP(GLS, ACO), KP, CVRP(ACO), MKP(ACO), BPP(ACO), CAF Templates & Text TSP(GLS, ACO, POMO, LEHD), CVRP(ACO, POMO, LEHD), MKP(ACO), BPP(ACO), DPP(GA) / TSP(GLS), FSSP(GLS) / PA(SA), SymmW(SA), SkewW(SA), BTD(GA), EPA(SA), FR(DFS) Template & Text Some open mathematical construction problems Callable funcs & Text BPP, TSP(GLS), CAF, BBOB, MIS, CVRP, TSP, PMSP

integrates RL techniques to mitigate this [12], [33], [34]. However, all these LHHs are only good at designing simple functions, suffering from suboptimal framework designs. AlphaEvolve [9] attempts to mitigate this issue by reducing token consumption through LLM-generated modification commands, yet it merely addresses the symptom (token pressure) rather than the root cause (structural limitations). Moreover, the conventional evaluation method for LHHs are casual, with reported improvements primarily stem from modifying trivial functions (sometimes even just a simple function returning the maximum of an array [7] or a standardized matrix [12] can get the best effects) within suboptimal algorithmic frameworks (often far from SOTA), artificially inflating their perceived capability while offering limited real-world applicability and scientific value. Bi-level Optimization. Bi-level Optimization (BLO) is a branch of mathematical programming [15] widely used in the Neural Architecture Search (NAS) field, implemented via evolution-based [35], gradient-based methods [36], [37]. Nevertheless, BLO applications in HH remain limited due to non-differentiability and large search spaces. Traditional HHs can only be considered as upper-level optimization frameworks that search for effective combinations of metaheuristic and parameter configurations to solve a given optimization problem [1], [38]. Even in works related to BLO, the outer layer is typically confined to hyperparameter tuning [39]. III. BEAM: B I - LEVEL M EMORY- ADAPTIVE A LGORITHMIC E VOLUTION As illustrated in Fig. 2, BEAM is composed of: 1) a core bi-layer structure inspired by modular programming [17], [44] (see Section III-A and Section III-B); 2) an external optimiza-

Single function Single function Single function Single function Entire algorithm Single function Single function

Single function Single function Entire algorithm Entire algorithm & Hybrid algorithm

tion mechanism called Adaptive Memory (see Section III-C). We also streamline LHHs and compare them detailedly in Table I. All the prompts used are provided in Appendix H. In this paper, instead of treating a complete algorithm as a single entity, we decompose it into a structure and function components to address the challenge of designing complete heuristics. Let I denote a heuristic individual (a complete algorithm), consisting of an algorithm structure S(I) and a set of functions F(I) = {f1 , . . . , fN }. We measure the overall quality Q(I) as the average performance (e.g., solution optimality gap or objective value) of individual I evaluated on a validation set of problem instances. This overall quality decomposes into two components: the structure quality Qs (S(I)) represents the inherent effectiveness of the overall algorithmic framework, not including the function implementations, and the function quality Qfi (fi | S(I)) measures how well each specific function fi implements its role within the given structure (e.g., a neighborhood evaluation function in local search). With the abovementioned definition, a bi-level formulation is a must since the quality of a heuristic cannot be determined until all its componentsP are implemented and executed N together: Q(I) = Qs (S(I)) + i=1 Qfi (fi | S(I)), where N is the number of functions required by the structure S(I). Thus, we formulate the bi-level optimization problem as follows: min Q(α, w∗ (α)),

(1)

α

s.t.

w∗ (α) = arg min Q(α, w), w

(2)

where Eq. (1) optimizes the structure, and Eq. (2) optimizes function realizations for a given structure. The upper-level

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

variable α is a symbolic representation (encoded as prompts and code templates for LLM) of the algorithm structure, corresponding to the Exterior Layer. The lower-level variable w represents the specific function implementations for a given structure α, corresponding to the Interior Layer, and w∗ (α) is the best realization conditioned on α. Note that both α and w are discrete symbolic objects in our LLM-based framework, though we use continuous notation for consistency with bilevel optimization literature. To balance exploration and exploitation, we optimize the Exterior Layer using a Genetic Algorithm (GA) and the Interior Layer using Monte Carlo Tree Search (MCTS) to solve the bi-level problem. The GA evolves the algorithm structures, while MCTS efficiently searches for high-quality function implementations within each structure. Details of the approach are described in the following sections. We use min because Q(·) represents solution quality measured as optimality gap or cost (lower is better); for maximization problems, objective values are negated before evaluation.

4

Algorithm 1 Genetic Evolutionary Process for Heuristic Individual Structure in BEAM (Exterior Layer). 1: t ← 0 2: P (0) ← LLM(task-prompt) ▷ Initialize population via LLM 3: P (0) ← Select(Education(P (0) ))) ▷ Educate and select 4: while t < T do 5: Sort P (t) by quality Q(I) in descending order 6: for i = 0 to len(P (t) ) − 2 do 7: if Bernoulli(pc ) = 1 then (t) (t) ▷ Cross consecutive individuals 8: Iˆcross ← Crossover(Ii , Ii+1 ) (t) (t) 9: Pcross ← Pcross ∪ {Iˆcross } 10: end if 11: end for 12: for i = 0 to len(P (t) ) − 1 do 13: if Bernoulli(pm ) = 1 then (t) ▷ Elitist mutation using best individual 14: Iˆmut ← Mutation(Ii ) (t) (t) 15: Pmut ← Pmut ∪ {Iˆmut } 16: end if 17: end for (t) (t+1) (t) 18: P ← Select(Education(P (t) ∪ Pcross ∪ Pmut )) 19: t←t+1 20: if t mod am_interval = 0 then ▷ Trigger AM every am_interval generations

21: AM ← UpdateAM(AM, P (t) ) 22: end if 23: end while 24: return Best individual from P (T )

▷ Algorithm 3

B. Interior Layer A. Exterior Layer We employ Genetic Evolution [30], [45] to evolve heuristic structures, following recent LHH literature [7], [11]. A popu(t) (t) (t) lation at generation t is P (t) = {I1 , I2 , . . . , In }, where (t) each Ii is a complete heuristic individual with structure (t) (t) S(Ii ) and functions F(Ii ). The population is updated through Algorithm 1. Individuals are sorted by quality Q(I) and processed in this quality-ranked order during crossover and mutation operations. Below we describe the evolutionary operators. Population initialization. This sector initializes P (0) by prompting the LLM with task descriptions, function signatures, requirements, HeuBase and KnoBase (see Section IV-B). Education & Selection. Before selection, BEAM first educates the population through the Education operation, which sends each structure to the Interior Layer (Section III-B) to realize its functions via MCTS, completing the individual. Then, the individuals are sorted according to Q(I) and if the population reaches max_pop_size, the worst-performing individuals will be eliminated. Crossover. For the implementation of Crossover(·, ·) in the prompt level, we simplify ReEvo’s approach by: 1) eliminating the resource-intensive reflection process [41], instead directly comparing solutions, and 2) restricting crossover to algorithm structure only, excluding functions to ensure meaningful comparisons. Mutation. We also followed ReEvo and used Elitist mutation [46] on the prompting strategy for Mutation(·), which requires LLM to learn from the best heuristic individual and redesign the current one. This operator is also performed on algorithm structure.

The Interior Layer implements the Education Operation, which completes and evaluates structures proposed by the Exterior Layer. Given a partial structure S(I) with placeholders for required functions F(I) = {f1 , . . . , fN }, Education realizes these functions, repairs generated code, and calibrates hyperparameters to produce a runnable individual I. The general process is presented in Algorithm 2. Monte-Carlo Tree Search (MCTS). Generally, the MCTS method has the valuation function1 : V (st ) ← V (st ) + α[rt+1 + γV (st+1 ) − V (st )]. Specifically, considering that the number of functions to design is finite and fixed for a certain individual, we use the recursion function: V (ft ) = max[r(ft ) + Vft (ft+1 )], ft

(3)

where V (ft ) =Qs (I) + Qft (I) + · · · + QfN (I), r(ft ) =Qft (I),

(4)

Vft (ft+1 ) =Qs (I) + Qft+1 (I) + · · · + QfN (I) with the t-th function being ft . Note that here we have Qs (I) ∈ α for all individuals I. Therefore, in order to choose the best function set of a heuristic individual, we need to maximize V (f1 , α), and w(α) = max V (F1 , α). In practice, for each function in the given structure, we try several different realization of the function, and then fill in all the functions that remains unrealized. After the evaluation of the entire structure, we select the best one’s function realization. This process will loop until all the functions are properly realized. 1 An alternative is One-Shot method, which simultaneously fills in all the functions. The choice between methods depends on both time constraints and the specific problem requirements.

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

Fixing. We add a fixing process [47] after the function fill-in process following LLaMEA [8] to address code errors. BEAM specifically handles: 1) compile/runtime errors, and 2) constraint violations, ensuring full heuristic exploitation. Calibration. Following prior work [41], [43], we add a calibration (hyperparameter tuning) feature: we require LLM to give a hyperparameter test range and then utilize the traditional technique CMA-ES [48] for calibration.

5

AM also maintains a fixed capacity Cmax and evicts lowutility entries using: ¯ ), U ∗ (f ) = λS(f ) + (1 − λ)∆(f

(6)

Algorithm 2 Monte-Carlo Tree Search with Fixing and Calibration for Function Realization in BEAM (Interior Layer).

¯ ) is the exponential moving average of recent where ∆(f fitness improvements. When capacity is exceeded, the lowestutility functions are removed, and rarely-used entries are pruned periodically. The algorithm is shown in Algorithm 3. Key notation: C is the candidate set from top elites, S(f ) is the selection score, and U ∗ (f ) is the long-term utility used for eviction decisions.

1: for t = 1 to N do ▷ N = total number of functions 2: FuncList ← LLM(I, Fill-t-prompt) ▷ Generate m candidates for ft ∗ 3: I ← LLM(I, FuncList0 , Fill-all-prompt) ▷ Complete with

Algorithm 3 Adaptive Memory insertion and maintenance.

ft = FuncList0 4: I ∗ ← Fix(I ∗ ) 5: best idx ← 0 6: for j = 1 to m − 1 do ▷ Try remaining candidates 7: I temp ← LLM(I, FuncListj , Fill-all-prompt) temp temp 8: I ← Fix(I ) 9: if Q(I temp ) > Q(I ∗ ) then ▷ Higher quality is better ∗ temp 10: I ←I ▷ Update best individual 11: best idx ← j 12: end if 13: end for ▷ Fix ft to best candidate 14: I ← I with ft ← FuncListbest idx 15: end for 16: I ← Calibration(I) ▷ Tune hyperparameters 17: return I

C. Adaptive Memory We also introduce a mechanism called Adaptive Memory (AM) to facilitate complex code generation. AM happens when the evolutionary process needs to reset population [49]. We present a new way beyond simply retaining elite individuals while injecting random new members [50]. The motivation behind AM is to make previously generated, high-quality functions directly reusable by the LLM. In LLMbased heuristic design, repeatedly generating long code blocks for every new candidate is costly and noisy. Instead of reducing output size by ‘patching’ existing code like AlphaEvolve [9], BEAM lets the LLM recall and call stored functions via importing them, greatly shrinking the amount of code it must emit on each generation. This approach also promotes diversity: by exposing a pool of reusable components, AM encourages the LLM to combine them in new ways just like how innovation often happens in scientific research. In practice, AM selects functions generated through MCTS that appear in elite solutions at the end of every am_interval generations (see Algorithm 1), and updates the pool iteratively by inserting strong new functions and retiring outdated ones. AM only provides the LLM with the stored functions’ names and purpose statements, allowing the LLM to directly import them. For each candidate function f ∈ C, we compute a composite score: S(f ) = α1 F̃fit (f ) + α2 F̃nov (f ) + α3 F̃use (f ) − α4 F̃age (f ), (5) with novelty measured by Fnov = 1 − maxg∈AM sim(f, g). If f is very similar to an existing memory entry g ∗ (similarity > τ ), it only replaces g ∗ when its score exceeds S(g ∗ ) + ∆th .

1: C ← extract functions from top-E elite individuals 2: for each f ∈ C do 3: Compute S(f ) using Eq. (5) 4: Find most similar function: g ∗ = arg maxg∈AM sim(f, g) 5: if sim(f, g ∗ ) > τ then ▷ High similarity detected 6: if S(f ) > S(g ∗ ) + ∆th then ∗ 7: Replace g with f in AM 8: else 9: Discard f ▷ Not sufficiently better 10: end if 11: else 12: Insert f into AM ▷ Novel function 13: end if 14: end for 15: while |AM| > Cmax do ▷ Capacity overflow 16: Evict function with smallest U ∗ (f ) (Eq. (6)) 17: end while 18: Remove functions unused for Tidle generations with U ∗ (f ) < ε

Note that this mechanism is uniquely enabled by our framework: traditional LHHs cannot support it due to their singlelayer structure. In all, AM allows LLM to retain high-performing functions from previous generations without realizing them again while encouraging new combinations. IV. K NOWLEDGE AUGMENTATION P IPELINE FOR LHH E VALUATION This section presents our Knowledge Augmentation (KA) pipeline for evaluating LHHs on complex code generation and complete solver design. We first discuss limitations in existing LHH evaluation practices (Section IV-A), then introduce our KA pipeline and its intended evaluation focus (Section IV-B). A. Limitations of Existing LHH Evaluation Current LHH benchmarks in the optimization field suffer from a fundamental mismatch with AHD requirements: 1) Most of them focus on designing small functions within predefined algorithms for problems like CO [11], [12]. They strongly rely on the human-designed external solver frameworks to achieve good results. 2) Most of them evaluated LHH-designed individual functions within suboptimal CO solvers [7], [41]. Reported improvements often stem from modifying trivial functions—in some cases, even a simple function computing array maxima [7] or normalizing matrices [12] is sufficient to achieve SOTA results. Consequently, these benchmarks offer limited practical utility. Moreover, the metric is also arbitrarily defined. Details can be found in Table I.

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

6

prior expert knowledge, which is then fed back into the LLM context. This makes the heuristic design process more informed and reduces reliance on the LLM’s latent parameter memory. V. E XPERIMENTS

Fig. 3: Proposed KA pipeline. For each task, LLMs first summarize problem-specific tags and then use them to construct KnoBase and the LLM-retrieved part of HeuBase.

B. KA-Guided Evaluation The need for KA comes from two observations about LLMs in heuristic design: • LLMs struggle to create very complex heuristics from scratch; they often either get stuck at initial solutions or generate constraint-violating outputs. • Yet LLMs are strong at repurposing components and at designing meta-frameworks such as outer architectures, parameter schedules, and heuristic coordination. Therefore, the core idea is to let LHHs design complete solvers: it evaluates how well an LHH composes reusable components, integrates domain knowledge, and produces working algorithms without relying on externally fixed solver frameworks. With this idea, Fig. 3 shows our KA pipeline, where LLMs build two structured databases: a HeuBase of callable functions and a text-based KnoBase of retrieved prior knowledge. a) HeuBase: HeuBase is a reusable heuristic component repository that LLMs can directly call, unlike AlphaEvolve’s template-based Program Database [9], [29]. We provide only function names and brief descriptions, similar to AM, and show that this supports novel combinations while preserving diversity (See B). HeuBase has two parts: LLM-Retrieved components, which are pip-installable libraries identified by LLM researchers from task tags and added via requirements.txt; and Pre-Constructed components, which are handcrafted heuristic routines (e.g. optimized 2-opt variants [51]) unavailable from pip packages. The database is intentionally compact but designed to grow through community contribution. Notably, while our KA pipeline focuses on designing complete solvers by composing reusable components, it establishes a complementary relationship with existing single-function LHHs. These LHHs excel at iteratively refining individual functions within fixed algorithmic frameworks, and their optimized single functions can be seamlessly integrated into HeuBase, enriching the repository and enhancing the overall design capabilities of complete-solver LHHs like ours. b) KnoBase: KnoBase is a text-based knowledge repository constructed from task-specific search tags produced by the LLM. LLM researchers use these tags to gather and summarize

In our experiments, we fixed the budgets for LHHs during the evolving stage. Generally, we strictly adhere to their original hyperparameter configurations, maintaining identical proportional relationships while only scaling magnitudes to ensure equivalent evaluation budgets. Besides, all experiments use consistent LLM models and temperature settings. However, due to the various difficulties of different problems, the budgets for different problems are different, as shown in Table II. Hardware and hyperparameter details are shown in Appendix E0g. For a given problem, we cover the test instances across all sizes with the same algorithm. The BEAM-generated algorithms are provided in Appendix J0d. A. On Traditional Single Function Evaluation In this section, we compare BEAM with other LHHs and expert-designed heuristics using traditional evaluation settings, averaging the performance of their best-generated heuristics over 5 trials. We use the three most-tested single function design tasks: 1) Design penalty heuristics in the Guided Local Search framework for the Traveling Salesman Problem (TSP). 2) Design priority functions for the Bin Packing Problem (BPP). 3) Design Cost-aware Acquisition Functions (CAF) for Bayesian Optimization (BO). For CO, we compare BEAM with ReEvo, EoH and MCTS-AHD. For CAF, we compare our results with more LHHs and expert-designed EI-cool [52] (EI [53] + EIpu [54]). Main Results. As shown in Table III, BEAM demonstrates strong overall performance while showing slight tendencies of overfitting in TSP and BPP. Note that EoH is worse than their published results [11] since we control the budget for running EoH. While their paper shows a final heuristic with 0.6% gap on Weibull5k, we cannot reproduce the result even if we triple the budget (> 2% gap). On CAF benchmarks (Table IV), BEAM outperforms AlphaEvolve in most datasets within same evolve budget. Other LHHs in the CAF experiment aren’t reimplemented and we simply test their best heuristic in their repository, so the budget is unknown for those LHHs. While our framework isn’t designed for single-function generation tasks - and consequently may introduce unnecessary complexity for such tasks - it nonetheless delivers competitive results. B. On Proposed KA-guided Evaluation We conduct experiments using our KA-guided evaluation pipeline. For CO problems, we implement runtime control by requiring LLM-generated algorithms to include a time-checking mechanism. This is implemented via Python’s time.time() function with a timeout parameter passed to the generated code. Runtime budgets for different problem sizes are listed in the table captions. In this section, we report

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

7

TABLE II: Detailed Budget Settings for Comparison with LHHs (Evolving Stage). Section

Problem

LLM Model

Budget Type

Budget

Attempts

Metric

Section V-A Section V-A Section V-A

TSP BPP CAF

DeepSeek-V3 DeepSeek-V3 DeepSeek-V3

Time Time Time

40min 60min 100min

5 5 5

Best & Average Best & Average Best

Section V-B Section V-B Section V-B Section V-B Section V-B

MIS (w/ RLSA) MIS (w/ KaHIP & ARW) CVRP TSP BBOB

DeepSeek-V3 DeepSeek-V3 DeepSeek-R1 DeepSeek-V3 DeepSeek-V3

Time Time Time Time Time

4h 4h 5h 1h 18min

3 3 3 3 3

Best∗ Best Best Best Best

Section V-C CVRP DeepSeek-R1 Token 55w 15 Best & Average Section V-C TSP DeepSeek-V3 Token 6w 15 Best & Average ∗: In Section V-B, since we are comparing LHHs with SOTA solvers, we only report the best results for comparison.

TABLE III: Results on Traditional CO Benchmark (TSP & BPP Single Function Design). Methods

TSP-100 MIN↓ AVG↓

ReEvo [7] EoH [11] MCTS-AHD [33] BEAM (ours)

0.01% 8.39e-3% 0.02% 2.63e-3%

0.03% 0.01% 0.04% 6.77e-3%

TSP-500 MIN↓ AVG↓

Weibull5k MIN↓ AVG↓

Weibull10k MIN↓ AVG↓

Weibull100k MIN↓ AVG↓

1.00% 0.85% 0.99% 0.88%

2.83% 3.13% 4.25% 2.58%

2.71% 2.90% 4.06% 1.99%

2.38% 2.75% 3.88% 1.82%

1.07% 0.96% 1.12% 0.95%

3.36% 3.18% 4.25% 3.26%

3.33% 3.02% 4.06% 3.04%

3.07% 2.87% 3.88% 2.86%

TABLE IV: Results on Traditional CAF Benchmark. ∗ Best heuristic from their repository; MH: Manually-designed Heuristic; TH: Trigonometric-Hump, Styblinski: Styblinski-Tang, HM: Hartmann. Methods

Type

Griewank↓ C=12 C=120

Rosenbrock↓ C=12 C=120

Levy↓ C=12 C=120

C=12

TH↓ C=120

7.14 4.75 19.81 15.93 9.81

0.43 0.11 0.06 0.15 0.26

1.78 0.24 0.08 1.89 1.68

1.05e-3 1.67e-3 3.34e-3 3.06e-3 2.70e-4

EI-cool EvoCAF∗ MCTS-AHD∗ AlphaEvolve BEAM (Ours)

MH LHH LHH LHH LHH

0.78 1.06 0.56 0.89 0.49

0.18 0.13 0.22 0.19 0.13

Methods

Type

HM3D↓ C=12 C=120

EI-cool EvoCAF∗ MCTS-AHD∗ AlphaEvolve BEAM (Ours)

MH LHH LHH LHH LHH

1.36e-2 5.93e-2 1.39e-2 3.76e-2 1.27e-2

9.27e-5 2.51e-3 2.06e-3 3.45e-3 3.95e-4

1.85 0.05 0.48 0.14 0.41

9.14e-4 1.77e-3 2.67e-3 2.18e-3 6.99e-4

Styblinski↓ C=12 C=120 11.88 1.93 0.69 14.82 1.25

3.92e-3 0.02 7.10e-3 0.02 1.54e-3

Powell↓ C=12 C=120

Shekel↓ C=12 C=120

HM6D↓ C=12 C=120

Cosine8↓ C=12 C=120

144.62 70.36 15.78 82.51 39.62

8.83 9.12 7.46 8.90 8.79

1.25 1.84 1.46 1.03 1.00

1.18 1.57 0.84 0.79 0.97

6.91 0.04 0.18 6.32 2.29

7.25 0.47 0.15 7.82 5.61

0.06 0.01 0.45 0.10 0.11

0.38 0.03 0.06 0.42 0.39

the best results after three trials. The results are shown in Figure 11, Table V and Table VI.

Main Results. For CO problems, BEAM surpasses KaMIS without reduction operations, achieves results close to HGS, and outperforms existing LHHs, demonstrating its superiority in complex code generation. Note that in TSP, since the ACO in ReEvo’s repository is a general framework without 2-opt [51], a robust ACO framework integrated with 2-opt can easily surpass EACO-EDM. BEAM and EoH both implement 2-opt. However, EoH fails to consistently outperform EACOEDM across all datasets, suggesting its suboptimal ACO design. For Continuous Optimization problems, results on BBOB show that BEAM can also achieves near-SOTA performance in continuous domains. Further insights derived from BEAMdesigned solvers can be found in Appendix L.

Fig. 4: Best individual distribution. Label ‘Gap’ is the performance of the code generated by the models. It shows that the code generated by BEAM has the best average performance and the smallest variance among the three models. C. More Comparative Results a) Best Individual Distribution.: We execute each of the three LHHs 15 times on CVRP (a comparatively complex

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

8

TABLE V: Performance comparison on MIS, CVRP, and TSP. Methods

Type

RB-200-300 (t=5s) T∗ OBJ↑ GAP↓

RB-800-1200 (t=60s) T OBJ↑ GAP↓

T

SATLIB (t=60s) OBJ↑ GAP↓

KaMIS (ReduMIS, 60s)

GA

15k

20.09

0.00%

15k

43.00

0.00%

15k

425.95

0.00%

RLSA ReEvo w/ RLSA EoH w/ RLSA MCTS-AHD w/ RLSA BEAM w/ RLSA

SA LHH LHH LHH LHH

9k 75 75 75 75

19.92 19.99 20.05 20.01 20.05

0.85% 0.50% 0.18% 0.39% 0.19%

60k 150 150 150 150

39.79 40.79 41.21 41.13 41.65

7.47% 5.03% 4.04% 4.35% 3.05%

60k 150 150 150 150

411.81 423.61 424.20 423.96 424.24

3.32% 0.55% 0.41% 0.47% 0.40%

ARW LS 500k 20.09 0.00% 2m 42.68 0.75% 18m 425.51 KaMIS (EvoMIS) GA 15k 20.09 0.01% 15k 42.97 0.06% 15k 425.95 EoH w/ KaHIP&ARW LHH 15k 20.09 0.00% 15k 43.01 -0.03% 15k 425.95 MCTS-AHD w/ KaHIP&ARW LHH 15k 20.09 0.00% 15k 42.97 0.08% 15k 425.91 AlphaEvolve w/ KaHIP&ARW LHH 15k 20.09 0.01% 15k 42.92 0.19% 15k 425.69 BEAM w/ KaHIP&ARW LHH 15k 20.09 0.00% 15k 43.03 -0.06% 15k 425.95 *: We control T (the iterations of local search algorithms) to ensure a similar runtime for fair comparison. Methods

Type

CVRP-100 (t=20s) OBJ↓ GAP↓

CVRP-200 (t=60s) OBJ↓ GAP↓

CVRP-500 (t=300s) OBJ↓ GAP↓

HGS [55]

GA

15.56

19.63

37.15

0.00%

0.00%

0.10% -1.44e-3% -1.45e-3% 0.01% 0.06% -1.93e-5%

0.00%

Split [56] & LS∗ LS 15.65 0.56% 20.00 1.89% 38.56 3.80% ReEvo w/ Split & LS LHH 15.62 0.37% 19.79 0.79% 37.71 1.52% MCTS-AHD w/ Split & LS LHH 15.58 0.13% 19.74 0.56% 37.63 1.29% EoH w/ Split & LS LHH 15.59 0.22% 19.74 0.57% 37.55 1.09% BEAM w/ Split & LS LHH 15.57 0.09% 19.70 0.38% 37.47 0.86% *: We perform the two algorithms on random permutations. The runtime is controlled to match its counterparts. Methods

Type

TSP-50 (t=5s) OBJ↓ GAP↓

TSP-100 (t=15s) OBJ↓ GAP↓

TSP-500 (t=40s) OBJ↓ GAP↓

EACO-EDM [7]

ACO

5.73

0.00%

8.13

0.00%

19.80

0.00%

EoH [11] w/ EDM [7] BEAM w/ EDM [7]

LHH LHH

5.76 5.73

0.52% -0.10%

8.13 7.90

-3.7e-4% -2.83%

18.11 17.69

-8.53% -10.66%

TABLE VI: Different LHHs on BBOB evaluation. ∗ From LLaMEA’s repository [57] (we test the ERADS function which is claimed to be its best design).

Methods

Rastrigin GAP↓

Rosenbrock GAP↓

Sphere GAP↓

Ackley GAP↓

Griewank GAP↓

Average GAP↓

LLaMEA [8]∗ ReEvo [7] EoH [11] LlaMEA-HPO BEAM

0.995 0.002 10.519 1.512 0.026

0.000 4.785 14.804 0.419 0.000

0.000 0.000 0.543 5.5e-10 0.000

4.4e-16 1.1e-5 0.714 6.9e-5 0.000

0.007 1e-6 3.5799 0.001 0.007

0.201 0.957 6.032 0.386 0.007

TABLE VII: Ablation study on Adaptive memory, education, and model-generalization (TSP is tested on TSP-500, MIS is tested on RB 800-1200, CVRP is tested on CVRP-500, CAF is tested on Ackley & Rastrigin.) Adaptive Memory

TSP

CVRP

CAF

BEAM BE

-9.55% -8.12%

0.86% 0.89%

3.46% 4.41%

Education Method

MIS

CVRP

CAF

One-Shot MCTS

3.63% 3.05%

1.07% 0.86%

5.12% 8.17%

Model Generalization

TSP-50

TSP-100

TSP-500

Deepseek-V3 GPT 3.5 turbo GPT 4o mini

0.00% 0.38% 0.19%

0.00% 2.64% 0.13%

0.00% 6.60% 0.12%

task) using the same evaluation dataset and record their best fitness values from each run. From Fig. 4, we can see that the average peak performance of BEAM far exceeds that of ReEvo

and EoH. Beside, EoH shows the poorest stability confirming HSEvo’s findings [41], and BEAM achieves exceptional stability. b) Evolution Curve: We analyze the median-performing evolutionary processes from 15 runs, tracking how their fitness values scale with token counts. Figures 7 and 8 show the performance gap (y-axis, lower is better) versus cumulative token consumption (x-axis) for TSP and CVRP respectively. Each curve represents one LHH framework’s evolutionary trajectory, where points indicate when new individuals are generated and evaluated. Among them, BEAM has the greatest improving ability. However, due to its bi-layer structure, it needs the largest number of tokens to generate its first heuristic individual. In CVRP (a more complex task), the initial vacancy is because most of the initial codes suffer from execution errors. c) Differences of Generated Heuristic: As illustrated in Fig. 5, BEAM consistently produces the longest and most complex heuristics (same requirement prompts). The heuristics

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

Fig. 5: Heuristic Length Distribution. BEAM consistently produces the longest and most complex heuristics, which is one reason for its superior performance.

9

Fig. 7: Gap-Token Tendency (TSP).

are also more robust compared to other LHHs, with detailed code comparison provided in Appendix K. D. Ablation Study a) On Adaptive Memory: Table VII shows that our proposed Adaptive Memory is effective. It also enhances evolution stability. Fig. 6 compares the iteration curves of BEAM and BE (BEAM without AM) during the evolution process. We observe that BEAM not only achieves higher peak performance but also exhibits better stability. To quantify it, Table VIII reports the mean and variance of the gaps at the final iteration over 5 independent runs. The results show that BEAM achieves both a higher mean and a lower variance, implying more robust and reliable convergence. Methods

AVG↓

VAR↓

BEAM BE

3.46 4.41

0.01 0.19

TABLE VIII: BEAM vs. BE: Average and Variance of Gaps.

Fig. 8: Gap-Token Tendency (CVRP).

c) On KA: The LLM-retrieved part of HeuBase is proved effective in BBOB, with CMA-ES being installed and utilized (see J0d). For KnoBase tests, we selected a Parallel Machine Scheduling Problem (PMSP) variant reformulated by ZeroBubble [58] used by DualPipe [59] in training DeepSeek-V3. Different from traditional 1F1B [60], this reformulated problem is fairly difficult with loads of constraints (see D). With KnoBase, BEAM incorporated the Sweep Line Algorithm and generated more constraint-satisfying initial solutions. d) Model Generalization: To evaluate the dependency on LLM size, we test smaller models on TSP, with results shown in Table VII. The result shows that despite slight increase in gaps, our BEAM still generates high-quality code with small models, showing that the demand and performance of LLM is not the decisive part of our BEAM. VI. C ONCLUSION AND F UTURE W ORK

Fig. 6: BEAM vs. BE: Best Gap - Iteration. b) On Individual Education Method: We test on two individual education methods and the results are shown in Table VII. We disable calibration for fairness. MCTS outperforms One-Shot except in CAF, where the objective is easy and the importance of structure outweighs functions.

In this paper, we propose BEAM, a Bi-layer structure that separates the heuristic design into two layers: the exterior layer for high-level algorithmic design and the interior layer for detailed function implementation. This structure, integrated with MCTS-based function selection and Adaptive Memory, enables the generation of high-quality, complex heuristics for both continuous and combinatorial optimization problems. We also introduce a Unified KA pipeline for more meaningful LHH evaluation. Experiments show that BEAM outperforms existing LHHs and SOTA solvers. Future work includes expanding BEAM to more complex domains and exploring more efficient ways for KA.

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

R EFERENCES [1] C. L. Camacho-Villalón, T. Stützle, and M. Dorigo, “Designing new metaheuristics: Manual versus automatic approaches,” Intelligent Computing, vol. 2, p. 0048, 2023. [2] J. H. Drake, A. Kheiri, E. Özcan, and E. K. Burke, “Recent advances in selection hyper-heuristics,” European Journal of Operational Research, vol. 285, no. 2, pp. 405–428, 2020. [3] N. Huynh and B. Lin, “Large language models for code generation: A comprehensive survey of challenges, techniques, evaluation, and applications,” 2025. [4] J. Jiang, F. Wang, J. Shen, S. Kim, and S. Kim, “A survey on large language models for code generation,” ArXiv, vol. abs/2406.00515, 2024. [5] P. Sahoo, A. K. Singh, S. Saha, V. Jain, S. Mondal, and A. Chadha, “A systematic survey of prompt engineering in large language models: Techniques and applications,” 2025. [6] J. Wei, X. Wang, D. Schuurmans, M. Bosma, B. Ichter, F. Xia, E. H. Chi, Q. V. Le, and D. Zhou, “Chain-of-thought prompting elicits reasoning in large language models,” in Proceedings of the 36th International Conference on Neural Information Processing Systems, ser. NIPS ’22. Red Hook, NY, USA: Curran Associates Inc., 2022. [7] H. Ye, J. Wang, Z. Cao, F. Berto, C. Hua, H. Kim, J. Park, and G. Song, “Reevo: Large language models as hyper-heuristics with reflective evolution,” in Advances in Neural Information Processing Systems, 2024. [8] N. v. Stein and T. Bäck, “Llamea: A large language model evolutionary algorithm for automatically generating metaheuristics,” IEEE Transactions on Evolutionary Computation, vol. 29, no. 2, pp. 331–345, 2025. [9] A. Novikov, N. Vũ, M. Eisenberger, E. Dupont, P.-S. Huang, A. Z. Wagner, S. Shirobokov, B. Kozlovskii, F. J. R. Ruiz, A. Mehrabian, M. P. Kumar, A. See, S. Chaudhuri, G. Holland, A. Davies, S. Nowozin, P. Kohli, and M. Balog, “AlphaEvolve: A coding agent for scientific and algorithmic discovery,” arXiv preprint arXiv:2506.13131, 2025. [10] F. Liu, R. Zhang, Z. Xie, R. Sun, K. Li, X. Lin, Z. Wang, Z. Lu, and Q. Zhang, “Llm4ad: A platform for algorithm design with large language model,” 2024. [11] F. Liu, X. Tong, M. Yuan, X. Lin, F. Luo, Z. Wang, Z. Lu, and Q. Zhang, “Evolution of heuristics: Towards efficient automatic algorithm design using large language model,” in International Conference on Machine Learning (ICML), 2024. [12] C. Mu, X. Zhang, and H. Wang, “Planning of heuristics: Strategic planning on large language models with monte carlo tree search for automating heuristic optimization,” 2025. [13] S. Lamm, P. Sanders, C. Schulz, D. Strash, and R. F. Werneck, “Finding near-optimal independent sets at scale,” 2015. [14] N. Pillay and R. Qu, Introduction to Hyper-Heuristics. Cham: Springer International Publishing, 2018, pp. 3–5. [15] B. Colson, P. Marcotte, and G. Savard, “An overview of bilevel optimization,” Annals OR, vol. 153, pp. 235–256, 06 2007. [16] S. Yao, D. Yu, J. Zhao, I. Shafran, T. L. Griffiths, Y. Cao, and K. Narasimhan, “Tree of thoughts: Deliberate problem solving with large language models,” ArXiv, vol. abs/2305.10601, 2023. [17] W. Zheng, S. P. Sharan, A. Jaiswal, K. Wang, Y. Xi, D. Xu, and Z. Wang, “Outline, then details: Syntactically guided coarse-to-fine code generation,” in International Conference on Machine Learning, 2023. [18] DeepSeek-AI, D. Guo, D. Yang, H. Zhang, J. Song, R. Zhang, R. Xu, Q. Zhu, S. Ma, P. Wang, X. Bi, X. Zhang, X. Yu, Y. Wu, Z. F. Wu, Z. Gou, Z. Shao, Z. Li, Z. Gao, A. Liu, B. Xue, B. Wang, B. Wu, B. Feng, C. Lu, C. Zhao, C. Deng, C. Zhang, C. Ruan, D. Dai, D. Chen, D. Ji, E. Li, F. Lin, F. Dai, F. Luo, G. Hao, G. Chen, G. Li, H. Zhang, H. Bao, H. Xu, H. Wang, H. Ding, H. Xin, H. Gao, H. Qu, H. Li, J. Guo, J. Li, J. Wang, J. Chen, J. Yuan, J. Qiu, J. Li, J. L. Cai, J. Ni, J. Liang, J. Chen, K. Dong, K. Hu, K. Gao, K. Guan, K. Huang, K. Yu, L. Wang, L. Zhang, L. Zhao, L. Wang, L. Zhang, L. Xu, L. Xia, M. Zhang, M. Zhang, M. Tang, M. Li, M. Wang, M. Li, N. Tian, P. Huang, P. Zhang, Q. Wang, Q. Chen, Q. Du, R. Ge, R. Zhang, R. Pan, R. Wang, R. J. Chen, R. L. Jin, R. Chen, S. Lu, S. Zhou, S. Chen, S. Ye, S. Wang, S. Yu, S. Zhou, S. Pan, S. S. Li, S. Zhou, S. Wu, S. Ye, T. Yun, T. Pei, T. Sun, T. Wang, W. Zeng, W. Zhao, W. Liu, W. Liang, W. Gao, W. Yu, W. Zhang, W. L. Xiao, W. An, X. Liu, X. Wang, X. Chen, X. Nie, X. Cheng, X. Liu, X. Xie, X. Liu, X. Yang, X. Li, X. Su, X. Lin, X. Q. Li, X. Jin, X. Shen, X. Chen, X. Sun, X. Wang, X. Song, X. Zhou, X. Wang, X. Shan, Y. K. Li, Y. Q. Wang, Y. X. Wei, Y. Zhang, Y. Xu, Y. Li, Y. Zhao, Y. Sun, Y. Wang, Y. Yu, Y. Zhang, Y. Shi, Y. Xiong, Y. He, Y. Piao, Y. Wang, Y. Tan, Y. Ma, Y. Liu, Y. Guo, Y. Ou, Y. Wang, Y. Gong, Y. Zou, Y. He, Y. Xiong, Y. Luo, Y. You, Y. Liu, Y. Zhou, Y. X. Zhu, Y. Xu, Y. Huang, Y. Li, Y. Zheng, Y. Zhu,

10

Y. Ma, Y. Tang, Y. Zha, Y. Yan, Z. Z. Ren, Z. Ren, Z. Sha, Z. Fu, Z. Xu, Z. Xie, Z. Zhang, Z. Hao, Z. Ma, Z. Yan, Z. Wu, Z. Gu, Z. Zhu, Z. Liu, Z. Li, Z. Xie, Z. Song, Z. Pan, Z. Huang, Z. Xu, Z. Zhang, and Z. Zhang, “Deepseek-r1: Incentivizing reasoning capability in llms via reinforcement learning,” 2025. [19] Z. Wang, Z. Zhu, Y. Han, Y. Lin, Z. Lin, R. Sun, and T. Ding, “Optibench: Benchmarking large language models in optimization modeling with equivalence-detection evaluation,” 2024. [20] X. Jiang, Y. Wu, Y. Wang, and Y. Zhang, “Bridging large language models and optimization: A unified framework for text-attributed combinatorial optimization,” 2024. [21] Z. Shi, M. Fang, and L. Chen, “Monte carlo planning with large language model for text-based game agents,” in The Thirteenth International Conference on Learning Representations, 2025. [22] A. AhmadiTeshnizi, W. Gao, and M. Udell, “Optimus: Scalable optimization modeling with (mi)lp solvers and large language models,” 2024. [23] C. Jiang, X. Shu, H. Qian, X. Lu, J. Zhou, A. Zhou, and Y. Yu, “Llmopt: Learning to define and solve general optimization problems from scratch,” 2025. [24] Z. Yuan, M. Liu, H. Wang, and B. Qin, “Ma-gts: A multi-agent framework for solving complex graph problems in real-world applications,” 2025. [25] X. Jiang, Y. Wu, C. Zhang, and Y. Zhang, “DRoc: Elevating large language models for complex vehicle routing via decomposed retrieval of constraints,” in The Thirteenth International Conference on Learning Representations, 2025. [26] Y. Gao, Y. Xiong, X. Gao, K. Jia, J. Pan, Y. Bi, Y. Dai, J. Sun, M. Wang, and H. Wang, “Retrieval-augmented generation for large language models: A survey,” 2024. [27] R. C. Staudemeyer and E. R. Morris, “Understanding LSTM - a tutorial into long short-term memory recurrent neural networks,” CoRR, vol. abs/1909.09586, 2019. [28] X. Wu, Y. Zhong, J. Wu, B. Jiang, and K. C. Tan, “Large language model-enhanced algorithm selection: Towards comprehensive algorithm representation,” in Proceedings of the Thirty-Third International Joint Conference on Artificial Intelligence, IJCAI-24, K. Larson, Ed. International Joint Conferences on Artificial Intelligence Organization, 8 2024, pp. 5235–5244, main Track. [29] B. Romera-Paredes, M. Barekatain, A. Novikov, M. Balog, M. P. Kumar, E. Dupont, F. J. R. Ruiz, J. Ellenberg, P. Wang, O. Fawzi et al., “Mathematical discoveries from program search with large language models,” Nature, vol. 625, no. 7995, pp. 468–475, 2024. [30] A. Eiben and J. Smith, Introduction To Evolutionary Computing. Springer, 01 2003, vol. 45. [31] Y. Zhou, A. I. Muresanu, Z. Han, K. Paster, S. Pitis, H. Chan, and J. Ba, “Large language models are human-level prompt engineers,” 2023. [32] R. Zhang, F. Liu, X. Lin, Z. Wang, Z. Lu, and Q. Zhang, “Understanding the importance of evolutionary search in automated heuristic design with large language models,” in Parallel Problem Solving from Nature – PPSN XVIII, M. Affenzeller, S. M. Winkler, A. V. Kononova, H. Trautmann, T. Tušar, P. Machado, and T. Bäck, Eds. Cham: Springer Nature Switzerland, 2024, pp. 185–202. [33] Z. Zheng, Z. Xie, Z. Wang, and B. Hooi, “Monte carlo tree search for comprehensive exploration in llm-based automatic heuristic design,” 2025. [34] A. Surina, A. Mansouri, L. Quaedvlieg, A. Seddas, M. Viazovska, E. Abbe, and C. Gulcehre, “Algorithm discovery with llms: Evolutionary search meets reinforcement learning,” 04 2025. [35] E. Real, A. Aggarwal, Y. Huang, and Q. Le, “Regularized evolution for image classifier architecture search,” Proceedings of the AAAI Conference on Artificial Intelligence, vol. 33, 02 2018. [36] H. Liu, K. Simonyan, and Y. Yang, “Darts: Differentiable architecture search,” 2019. [37] X. Wang, Z. Lian, J. Lin, C. Xue, and J. Yan, “Diy your easynas for vision: Convolution operation merging, map channel reducing, and search space to supernet conversion tooling,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 45, no. 11, pp. 13 974– 13 990, 2023. [38] V. Nannen and A. Eiben, “A method for parameter calibration and relevance estimation in evolutionary algorithms,” in Genetic and Evolutionary Computation Conference, vol. 1, 07 2006, pp. 183–190. [39] H. Park, H. Kim, H. Kim, J. Park, S. Choi, J. Kim, K. Son, H. Suh, T. Kim, J. Ahn, and J. Kim, “Versatile genetic algorithm-bayesian optimization(ga-bo) bi-level optimization for decoupling capacitor placement,” in IEEE 32nd Conference on Electrical Performance of Electronic Packaging and Systems (EPEPS), 10 2023, pp. 1–3.

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

[40] Y. Yao, F. Liu, J. Cheng, and Q. Zhang, “Evolve cost-aware acquisition functions using large language models,” in Parallel Problem Solving from Nature – PPSN XVIII, M. Affenzeller, S. M. Winkler, A. V. Kononova, H. Trautmann, T. Tušar, P. Machado, and T. Bäck, Eds. Cham: Springer Nature Switzerland, 2024, pp. 374–390. [41] P. V. T. Dat, L. Doan, and H. T. T. Binh, “Hsevo: Elevating automatic heuristic design with diversity-driven harmony search and genetic algorithm using llms,” in The 39th Annual AAAI Conference on Artificial Intelligence, 2025. [42] N. van Stein, D. Vermetten, and T. Bäck, “In-the-loop hyper-parameter optimization for llm-based automated design of heuristics,” ACM Transactions on Evolutionary Learning and Optimization, Apr. 2025. [43] C. D. Rosin, “Using code generation to solve open instances of combinatorial design problems,” 2025. [44] W. Sun, X. Sun, and X. Wang, “Using modular programming strategy to practice computer programming: A case study,” in ASEE Annual Conference and Exposition, 06 2012. [45] X. Wu, S. hao Wu, J. Wu, L. Feng, and K. C. Tan, “Evolutionary computation in the era of large language model: Survey and roadmap,” 2024. [46] K. Liagkouras and K. Metaxiotis, “An elitist polynomial mutation operator for improved performance of moeas in computer networks,” in 2013 22nd International Conference on Computer Communication and Networks (ICCCN), 2013, pp. 1–5. [47] R. Tian, Y. Ye, Y. Qin, X. Cong, Y. Lin, Y. Pan, Y. Wu, H. Hui, W. Liu, Z. Liu, and M. Sun, “Debugbench: Evaluating debugging capability of large language models,” 2024. [48] N. Hansen and A. Ostermeier, “Completely derandomized selfadaptation in evolution strategies,” Evolutionary Computation, vol. 9, pp. 159–195, 06 2001. [49] A. Fukunaga, “Restart scheduling for genetic algorithms,” Lecture Notes in Computer Science, 05 2002. [50] C. Li, Y. Liu, Y. Zhang, M. Xu, J. Xiao, and J. Zhou, “A novel multilevel population hybrid search evolution algorithm for constrained multiobjective optimization problems,” Journal of King Saud University Computer and Information Sciences, vol. 34, no. 10, Part B, pp. 9071– 9087, 2022. [51] J. K. Lenstra and E. Aarts, Local Search in Combinatorial Optimization. Princeton University Press, 2003. [52] E. H. Lee, V. Perrone, C. Archambeau, and M. W. Seeger, “Cost-aware bayesian optimization,” CoRR, vol. abs/2003.10870, 2020. [53] J. Mockus, “On bayesian methods for seeking the extremum,” in Optimization Techniques, 1974. [54] J. Snoek, H. Larochelle, and R. P. Adams, “Practical bayesian optimization of machine learning algorithms,” in Advances in Neural Information Processing Systems, F. Pereira, C. Burges, L. Bottou, and K. Weinberger, Eds., vol. 25. Curran Associates, Inc., 2012. [55] T. Vidal, T. G. Crainic, M. Gendreau, N. Lahrichi, and W. Rei, “A hybrid genetic algorithm for multidepot and periodic vehicle routing problems,” Operations Research, vol. 60, pp. 611–624, 06 2012. [56] C. Prins, “A simple and effective evolutionary algorithm for the vehicle routing problem,” Computers & Operations Research, vol. 31, no. 12, pp. 1985–2002, 2004. [57] N. van Stein and T. Bäck, “Llamea,” Sep. 2024, accessed: YYYY-MMDD. [58] P. Qi, X. Wan, G. Huang, and M. Lin, “Zero bubble pipeline parallelism,” 2023. [59] DeepSeek-AI, A. Liu, B. Feng, B. Xue, B. Wang, B. Wu, C. Lu, C. Zhao, C. Deng, C. Zhang, C. Ruan, D. Dai, D. Guo, D. Yang, D. Chen, D. Ji, E. Li, F. Lin, F. Dai, F. Luo, G. Hao, G. Chen, G. Li, H. Zhang, H. Bao, H. Xu, H. Wang, H. Zhang, H. Ding, H. Xin, H. Gao, H. Li, H. Qu, J. L. Cai, J. Liang, J. Guo, J. Ni, J. Li, J. Wang, J. Chen, J. Chen, J. Yuan, J. Qiu, J. Li, J. Song, K. Dong, K. Hu, K. Gao, K. Guan, K. Huang, K. Yu, L. Wang, L. Zhang, L. Xu, L. Xia, L. Zhao, L. Wang, L. Zhang, M. Li, M. Wang, M. Zhang, M. Zhang, M. Tang, M. Li, N. Tian, P. Huang, P. Wang, P. Zhang, Q. Wang, Q. Zhu, Q. Chen, Q. Du, R. J. Chen, R. L. Jin, R. Ge, R. Zhang, R. Pan, R. Wang, R. Xu, R. Zhang, R. Chen, S. S. Li, S. Lu, S. Zhou, S. Chen, S. Wu, S. Ye, S. Ye, S. Ma, S. Wang, S. Zhou, S. Yu, S. Zhou, S. Pan, T. Wang, T. Yun, T. Pei, T. Sun, W. L. Xiao, W. Zeng, W. Zhao, W. An, W. Liu, W. Liang, W. Gao, W. Yu, W. Zhang, X. Q. Li, X. Jin, X. Wang, X. Bi, X. Liu, X. Wang, X. Shen, X. Chen, X. Zhang, X. Chen, X. Nie, X. Sun, X. Wang, X. Cheng, X. Liu, X. Xie, X. Liu, X. Yu, X. Song, X. Shan, X. Zhou, X. Yang, X. Li, X. Su, X. Lin, Y. K. Li, Y. Q. Wang, Y. X. Wei, Y. X. Zhu, Y. Zhang, Y. Xu, Y. Xu, Y. Huang, Y. Li, Y. Zhao, Y. Sun, Y. Li, Y. Wang, Y. Yu, Y. Zheng, Y. Zhang, Y. Shi, Y. Xiong, Y. He, Y. Tang, Y. Piao, Y. Wang, Y. Tan, Y. Ma, Y. Liu, Y. Guo, Y. Wu, Y. Ou,

11

Y. Zhu, Y. Wang, Y. Gong, Y. Zou, Y. He, Y. Zha, Y. Xiong, Y. Ma, Y. Yan, Y. Luo, Y. You, Y. Liu, Y. Zhou, Z. F. Wu, Z. Z. Ren, Z. Ren, Z. Sha, Z. Fu, Z. Xu, Z. Huang, Z. Zhang, Z. Xie, Z. Zhang, Z. Hao, Z. Gou, Z. Ma, Z. Yan, Z. Shao, Z. Xu, Z. Wu, Z. Zhang, Z. Li, Z. Gu, Z. Zhu, Z. Liu, Z. Li, Z. Xie, Z. Song, Z. Gao, and Z. Pan, “Deepseek-v3 technical report,” 2025. [60] D. Narayanan, A. Harlap, A. Phanishayee, V. Seshadri, N. R. Devanur, G. R. Ganger, P. B. Gibbons, and M. Zaharia, “Pipedream: generalized pipeline parallelism for dnn training,” in Proceedings of the 27th ACM Symposium on Operating Systems Principles, ser. SOSP ’19. New York, NY, USA: Association for Computing Machinery, 2019, p. 1–15. [61] Z. Sun and Y. Yang, “DIFUSCO: Graph-based diffusion solvers for combinatorial optimization,” in Thirty-seventh Conference on Neural Information Processing Systems, 2023. [62] S. Feng and Y. Yang, “Regularized langevin dynamics for combinatorial optimization,” 2025. [63] S. Sanokowski, S. Hochreiter, and S. Lehner, “A diffusion model framework for unsupervised neural combinatorial optimization,” 2024. [64] J. Ma, W. Pan, Y. Li, and J. Yan, “Coexpander: Adaptive solution expansion for combinatorial optimization,” in International Conference on Machine Learning (ICML), 05 2025. [65] D. Drakulic, S. Michel, and J.-M. Andreoli, “Goal: A generalist combinatorial optimization agent learner,” 2025. [66] P. Sanders and C. Schulz, “Kahip v3.00 – karlsruhe high quality partitioning – user guide,” 2020. [67] M. Ngisomuddin and D. Satyananda, “Perturbation operator analysis on ils-rvnd algorithm to solve cvrp,” in THE 3RD INTERNATIONAL CONFERENCE ON MATHEMATICS AND SCIENCE EDUCATION (ICOMSE), vol. 2215, 04 2020, p. 070015.

A PPENDIX A. Additional Comparison Our usage of MCTS is fundamentally different from PoH and MCTS-AHD as detailed in Table X. B. Diversity Discussions Counterintuitively, the introduction of HeuBase enhances rather than diminishes the diversity of LLM-generated functions. This phenomenon can be attributed to two primary reasons: 1) diverse application of components (e.g. using KaHIP for either initialization or crossover), and 2) stochastic selection Table IX shows the selection frequency and drop rate in an extra MIS experiment with 50 samples, showcasing diversity. TABLE IX: Selection Frequency of HeuBase Functions. Function kaffpa node_separator arw arw_1iter deterministic_rounding

Selection Frequency 82% 56% 80% 86% 14%

C. Evaluation Overview Table XI documents all evaluation configurations. Note that the first three are chosen from some most-used LHH experiments. We exclude experiments like designing heuristic guide functions for ACO solving problems like TSP and CVRP [7], [33] since empirical results show that there’s no obvious performance difference between LHHs. We also replace the TSP w/ EDM task (described in Section V-B) with the complete TSP solver design since EDM is a comparatively trivial function and LHHs can also design good TSP solver without it.

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

12

TABLE X: Further Comparison. Method MCTS-AHD PoH BEAM (Ours)

MCTS Usage This work just substitutes EoH’s evolutionary operator with an MCTS-like strategy for algorithm search, still being a single-layer framework. The only difference between PoH and MCTS-AHD is that PoH introduced a reflection step similar to ReEvo in the MCTS process. In our work, MCTS is used as the education strategy for an algorithm structure in the interior layer of our bi-layer structure, searching different realization for the unrealized functions of the structure.

TABLE XI: Unified LHH Evaluation Settings with KA. Prob.

MC

Type

Design Type

Heuristic Type

Dataset Size

Allowed KA

TSP∗

Easy

GCO

Single Function

Not Restricted

3*128

/

BPP

Easy

CO

Single Function

Not Restricted

3*5

/

CAF

Easy

BO

Single Function

Not Restricted

10*2*5

/

MIS

Medium

GCO

Hybrid Algorithm

GA

3*500

HeuBase (RLSA)

MIS

Medium

GCO

Hybrid Algorithm

GA

3*500

HeuBase (KaHIP, ARW)

TSP

Medium

GCO

Entire Algorithm

ACO

3*128

/

BBOB

Medium

BBO

Entire Algorithm

Not Restricted

5

HeuBase KnoBase

CVRP

Hard

GCO

Hybrid Algorithm

Not Restricted

3*100

HeuBase (Split, LS)

PMSP

Hard

CO

Entire Algorithm

Not Restricted

4*3

HeuBase KnoBase

(LLM-Retrieved),

(LLM-Retrieved),

*: The target is to design a penalty heuristic within the Guided Local Search framework with perturbation_moves set to 30 and iter_limit set to 1200; GCO: Graph Combinatorial Optimization, BO: Bayesian Optimization, BBO: Black Box Optimization, ACO: Ant Colony Optimization; MC: Model Complexity

D. Formal Problem Descriptions a) Traveling Salesman Problem (TSP).: Given a complete graph G = (V, E) with |V| = N nodes and a symmetric cost matrix C ∈ RN ×N where Cij = Cji denotes the cost of traveling between nodes i and j, the objective is to find a Hamiltonian cycle τ = (i1 , i2 , . . . , iN , i1 ) that starts and ends at the same node, visits allP other nodes exactly once, and N −1 minimizes the total tour cost: k=1 Cik ik+1 + CiN i1 , with Cii = 0 for all i ∈ V. b) Online Bin Packing Problem (BPP): Given a sequence of items {x1 , x2 , . . . , xN } with sizes xi ∈ (0, 1] arriving one by one, the goal is to assign each item to a bin upon arrival without knowledge of future items. Each bin has a capacity of 1, and no bin may exceed this capacity. The objective is to minimize the total number of bins used to pack all items. c) Cost-aware Acquisition Functions (CAF) for Bayesian Optimization (BO): In Bayesian Optimization (BO), we aim to optimize an unknown function f : X → R with evaluation cost c(x) > 0 varying over x ∈ X . A CAF is defined as αCAF (x) := α(x) c(x) where: α(x) is a standard acquisition function such as Expected Improvement (EI), Upper Confidence Bound (UCB), or Probability of Improvement (PI); c(x) is the cost of evaluating f at point x. The next query point is then selected by solving: xt+1 = arg maxx∈X αCAF (x) = arg maxx∈X α(x) c(x) which prioritizes locations that offer the highest expected gain per unit cost. d) Maximum Independent Set (MIS): Given a unweighted graph G = (V, E), an independent set S ⊆ V is a subset of nodes such that no two nodes in S are adjacent. The goal is to maximize |S| s.t. ∀i, j ∈ S, (i, j) ∈ / E.

e) Capacitated Vehicle Routing Problem (CVRP): Given a graph G = (V, E), a depot node v0 ∈ V, a cost matrix C ∈ RN ×N , a demand vector d ∈ RN + , and a vehicle capacity Q > 0, the goal is to plan a set of routes R, each route r ∈ R starting and ending at the depot v0 , such that each customer node is visited exactly once and P the total demand on each route does not exceed Q, i.e., i∈r di ≤ Q. P The objective is P Cij . to minimize the total cost of all routes: min R r∈R (i,j)∈r

f) Black Box Optimization Benchmark (BBOB): Black-Box Optimization Benchmarking (BBOB) is COCO’s standard suite of 24 noiseless, single-objective test functions—provided in dimensions 2, 3, 5, 10, 20, and 40 with multiple randomized instances—to objectively compare black-box optimizers under a fixed function-evaluation budget. Performance is measured by the number of evaluations required to reach target accuracies, convergence curves, and success rates across functions of varying separability, conditioning, and multimodality. g) Parallel Machine Scheduling Problem (PMSP): This problem is reformulated by ZeroBubble [58]: Any pass in a pipeline can be uniquely identified by a triple (i, j, c), where i ∈ {1, 2, . . . , p} denotes the stage, j ∈ {1, 2, . . . , m} denotes the microbatch index, and c ∈ {F, B, W } represents the computation type (forward, backward, weight update). T(i,j,c) is the execution time of pass (i, j, c), and E(i,j,c) is its ending time. ∆M(i,j,c) denotes the memory change incurred by pass (i, j, c). For example, ∆M(·,·,F ) = MB indicates that the forward pass increases memory usage by MB . Similarly, the backward pass frees MB while requiring memory for weights MW , hence ∆M(·,·,B) = MW − MB , and the weight

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

update consumes ∆M(·,·,W ) = −MW . A binary indicator O(i,j,c)→(i′ ,j ′ ,c′ ) equals 1 if pass (i, j, c) is scheduled before pass (i′ , j ′ , c′ ), and 0 otherwise. The PMSP is then formulated as a Mixed Integer Linear Programming (MILP) problem: min

max E(i,m,W ) − E(i,1,F ) + T(i,1,F )

s.t.

E(i,j,F ) ≥ E(i−1,j,F ) + Tcomm + T(i,j,F )

O,E

i

E(i,j,B) ≥ E(i+1,j,B) + Tcomm + T(i,j,B) E(i,j,c) ≥ E(i,j,c′ ) + T(i,j,c) − O(i,j,c)→(i,j,c′ ) · ∞ X Mlimit ≥ ∆M(i,j ′ ,c′ ) + ∆M(i,j,c) O(i,j,c)→(i,j ′ ,c′ ) j,c

13

F. Hyperparameter Setting Common hyperparameter setting is presented in Table XIII. In the ablation study on education methods, we bypass KA and calibration, and Table XIV shows other settings. In the ablation study on AM, we refresh the population for BE (BEAM without AM) by keeping 2 elite algorithms and injecting 13 newly sampled ones at the same intervals as BEAM. In the ablation study on KA, all settings are the same except whether to include KA. In this section, TSP is tested on TSP-500, CVRP is tested on CVRP-500, MIS is tested on RB-800-1200, and CAF is tested on Ackley and Rastrigin (average over the two).

(7) TABLE XIII: Common Hyperparameter Setting for BEAM. TABLE XII: Evaluation Dataset for BEAM. Problem

Dataset

Instances

TSP (GLS) BPP CAF MIS (w/ RLSA) MIS (w/ KaHIP & ARW) CVRP TSP (w/ EDM) BBOB PMSP

TSP-200 Weibull5k Ackley & Rastrigin RB 200-300 RB 800-1200 CVRP-100 TSP-50 Ellipsoidal & Levy Randomly-Generated

10 10 5*2 25 10 20 30 2 5

Parameter

llm_temperature 0.7 (for Fixing) / 1.0 (others) crossover_rate 0.7 mutation_rate 0.3 mc_func_pop 3 max_func_num 4 mc_func_pop: The number of functions generated during each func i generation.

TABLE XIV: Hyperparamter Setting for Education Method Ablation Study. Methods

E. Dataset Details Note that here only provides the test dataset. The evaluation dataset during the LHH process isn’t restricted for this setting. The evaluation dataset we use is presented in Table XII. a) TSP: We follow DIFUSCO [61] to conduct experiments on TSP-50, TSP-100 and TSP-500. b) BPP: Following FunSearch [29], we used instances sampled from Weibull distributions. Specifically, we generated Weibull 5k, Weibull 10k, Weibull 100k. c) CAF: Following MCTS-AHD [33], all the dataset are synthetic instances with different landscapes and input dimensions, and we used Ackley and Rastrigin as the evaluation dataset during evolution, so these two aren’t included in the table below. We tested with sampling budgets of 12 and 120 to assess the generalizability of all algorithms. All tests take 5 trials. d) MIS: We use the Revised Model B (RB) graphs and SATLIB graphs, following [62], [63]. For RB graphs, we use RB 200-300 for small-scale and RB 800-1200 for large-scale. e) CVRP: We constructed CVRP-100, CVRP-200 and CVRP-500. Following COExpander [64] and GOAL [65], where the coordinates of the depot and clients were sampled from a uniform distribution over the unit square, consistent with the TSP setting. f) BBOB: We choose five functions in BBOB: Rastrigin, Rosenbrock, Sphere, Ackley and Griewank, and use the provided extrema points to evaluate. g) PMSP: All the data used for evaluation are real-world data directly taken from the ZeroBubble literature [58].

Value

iter

ips

mps

mft

One-Shot 8 20 5 3 MCTS 3 5 3 3 ips: init_pop_size, mps: max_pop_size,mft : max_fix_try

G. Hardware Details All experimental evaluations, including both the evolutionary optimization process and final performance assessments, are conducted on an Apple M3 CPU. However, algorithms incorporating RLSA are executed on an NVIDIA® GeForce RTX™ 4070 Ti SUPER GPU due to their PyTorch-based computational requirements and algorithms with KaHIP are run on an Intel(R) Xeon(R) Platinum 8558 96-Core Processor CPU since KaHIP doesn’t support ARM64 architecture. H. More Details on MCTS Fig. 9 further illustrates the MCTS process. Note that in practice, when prompting the LLM to generate multiple variants for a certain function, we provide its previous designs and tell LLM to “improve it” or “think about a different way to implement this” to encourage diversity (See Appendix H for details). Across multiple runs and different problems, MCTS granted an average performance gain of 46.6% to each individual. I. Common Prompts a) am.txt: This prompt is used to require the LLM to provide a name and a description for the given function, which will be later added to the Adaptive Memory. { function code}

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

14

Fig. 9: MCTS Example.

Please read the function code above carefully, and understand what this function achieves. Your job is to give a name of this function and then provide a short description for this function. Note that the description should include the overall description, the meaning of the arguments and the meaning of the output. Your output should follow this: ```python def [the function name] (...copy the provided arguments): ””” [the overall description] Args: [arg name]: [the meaning of the arguments] Out: [output argument name]: [the meaning of the output] ””” ``` b) ask pms interval.txt: This prompt is aimed to require LLM to provide the range of hyperparameters for Calibration part.

A heuristic will be given below, all the hyperparameters will be on the top. You should read the entire code carefully and decide on an interval for each of the hyperparameter. You should output a python dictionary only, and the dictionary maps from a string (name of the hyperparameter) to a tuple (start, end). The dictionary must be named pms dict. Format your dictionary as a Python code string: ”``p̀ython ... ```” c) ask pms system.txt: This prompt is the system prompt of ask pms interval. You are an expert in hyperparameter optimization and you give proper test range for each hyperparameter. Your suggestions on test range should be in the form of python dictionary. Your response outputs Python code and nothing else. Format your code as a Python code string: ”``p̀ython ... ```”.

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

d) crossover.txt: This prompt is the crossover prompt. the {} part will be fill in with the parent-structures. {exterior user generator} [Worse code] {worse code} [Better code] {better code} [Improved code] Please reflect on why the latter one perform better and write an improved structure according to your reflection. Enclose your code with a Python fenced code block. e) exterior user generator.txt: This prompt requires LLM to design the overall structure of the problem. Some requirements and restrictions are provided. Note that we call both AM and HeuBase are called ’Heuristic Database’ in the prompt to make a clearer understanding. Now you have to design a novel {alg type} solving {problem}. {problem description} You need to design the overall structure of your {alg type} as well as thinking detailedly about what you will do in each step and make sure each goal is achievable. Then you must write a piece of code, presenting your algorithm. Here are the requirements of your python code: 1. Put your main structure in the following function: {baseline} 2. You should look at the *heuristic database* first, and then think: - How can you deconstruct the big problem into small subproblems step by step? - What are the main goal of the subproblems that are needed to complete the heuristic designing? Your thoughts should guide you to complete the heuristic design in step 3. You should not think too much of how to realize the subproblems. 3. In this code you needn’t implement everything and utilize the idea of **modularization programming**. You can call external functions to represent or compose every subproblems. There are two cases of this external function: (The first case has a higher priority) I. The function already exists in the *heuristic database* (which I’ll give you later). You can suppose I’ve already implemented it and directly call it. You must make sure that the function 100% satisfies your need here. You’re encouraged to integrate the existing heuristics into your structure. **Note**: - I’ll import these heuristics when I test your code, so don’t define these heuristics yourself! II. The function doesn’t exist in the *heuristic database*. Remember to name these external functions ’func id’ where id is a **number**. It should

15

format as this: ```python def func {{id}}(...) -¿ ... : # Purpose: pass ``` **Note**: - You MUSTN’T implement these functions since I’ll let others implement it. - The purpose should be very clear. - Don’t bother calling these functions if the purpose is too easy! - There should be at least 1 func{{id}}, at most max func pop func{{id}}. The {{id}} must start from 1. 4. Put the definitions of all the hyperparameters on top of everything. Enclose your hyperparameter list with two ”#Hyperparameter#”. Your hyperparameters must be float or int. If a hyperparameter is int, add ”# int” right after the definition inline. 5. We only have {timeout} seconds to perform the algorithm, so you must set your code a timeoutsecond-clock. Include MAX TIME = {timeout} in your hyperparameter list. In the code, please make frequent check whether the time is up. 6. Let your code print out the best objectives after each iteration. {prior knowledge} f) fill 1func.txt: This prompt is used in MCTS, where LLM need to temperately fill in only one function to choose the best one to really fill into the structure. An algorithm solving {problem} will be given below, with some of the functions realized and others unrealized. {problem description} You are required to complete func {id}. You should follow the instructions given. You should output only python code as required. Your generated code MUST be different from the following code and you should either improve it or think out of box and explore a different way: {code before} **Critical Reminder:** - You MUST keep ALL existing code, comments, and hyperparameters EXACTLY as provided - Your response MUST contain the ENTIRE algorithm code - Only modify the specified func {id} implementations - Preserve ALL other code exactly as provided - Format output as: ```python [COMPLETE CODE] ``` prior knowledge g) fill allFunc.txt: This prompt is used in One-Shot method and MCTS (used after fill 1func to help decide which is the best function to actually fill into the structure). An algorithm solving **{problem}** will be given

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

below, with some of the functions realized and others unrealized. {problem description} You have to implement all the functions named func i() where i is a number. You should follow the instructions given in the code structure when you implement the functions and make sure your code serves the original purpose. You are encouraged to directly call the functions in the *heuristic database* that will be given below to serve your purpose. **Critical Reminder:** - You MUST keep ALL existing code, comments, and hyperparameters EXACTLY as provided - Your response MUST contain the ENTIRE algorithm code - Only modify the specified func i() implementations - Preserve ALL other code exactly as provided - Format output as: ```python [COMPLETE CODE] ``` {prior knowledge} h) fix.txt: This prompt is used to require LLM to fix the heuristic that reported error. The following code has the following error: {error msg}. Please fix it. **Important**: - You must **output the entire code** and you must **only** fix the error in the traceback message. Don’t fix anything besides the error presented by the error message. - You can’t change MAX TIME. - Don’t remove hyperparameters! i) fix system.txt: This prompt is the system prompt of fix. You are an expert in debugging and you output the entire code after debugging. Your response outputs Python code and nothing else. Format your code as a Python code string: ”```python ... ```”. j) func generation.txt: This prompt is the system prompt of all of the function generation process. You are an expert in heuristic design. Your task is to complete the functions in the given heuristic structure to meet the following requirements: 1. Your design must fully satisfy the requirements specified in the provided function templates. 2. Your design must maintain the same parameters as the given function templates. 3. Your design should work correctly within the context of the heuristic, which will be provided below. Your response must include the complete heuristic code with all necessary functions filled in. Your response outputs Python code and nothing else. Format your code as a Python code string: ”```python ... ```”.

16

k) heubase common.txt: This prompt is used to provide LLM with AM as well as the HeuBase. Below is the heuristic database. You can directly call any of them. (Note that ‘edge index‘ is sized (2, num edge) for unweighted graph and is sized (3, num edge) for weighted graph where the third dimension is the length of the edge; ‘ini sol‘ must be valid) Caution! Do not reimplement these functions—they are preloaded and conflicts will arise if duplicated. Simply call them by name with the required arguments. l) mutation.txt: This prompt is the mutation prompt. the {} part will be fill in with the relative-structures. {exterior user generator} [Now Structure] {now structure} [Elitist Code] {elitist structure} [Improved code] Please write a mutated structure based on the Now Structure. You should reflect on why the elitist code perform the best and take inspiration from it. Enclose your code with a Python fenced code block. m) pip search.txt: This prompt is used to require the LLM to search for the relative libraries related that may help construct the heuristic. We adopt Lepton AI2 for online search. We are solving a problem of problem name, problem description. You are required to search for some libraries for python that has a close relationship with this problem in topics or in details. Your output should follow this, act like a requirements.txt: ``` [library name 1] == [version number] [library name 2] == [version number] ``` n) prior knowledge.txt: This prompt is used to provide LLM with KnoBase. You may refer to these prior expert knowledge: {prior knowledge} o) problem description.txt: This prompt is used to provide LLM with problem description. The problem description is as follows: {problem description} p) system generator.txt: This prompt is the system prompt of the whole process. You are an expert in the domain of optimization heuristics. Your task is to design heuristics that can effectively solve optimization problems. Your response outputs Python code and nothing else. Format your code as a Python code string: ”```python ... ```”. 2 https://github.com/leptonai/search with lepton

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

J. Prompt Format for Specific Problems a) description.txt: This is different for every problem. Note that this is optional, because when facing a new problem, LLM needs to know the definition of the problem, while facing an old problem, there is no such need. For each problem, the description of the problem will be provided here. b) function signature.txt: This is different for every problem. For each problem, a function signature will be provided to specify the required input and output format. ```python def heuristic (...parameters...) -¿ ... ””” Args:... Returns:... ””” ``` c) heubase.txt: This is different for every problem. For each problem, this will provide LLM with Heubase and Adaptive Memory function’s conclusive summary, so that LLM can know the function’s usage without the need to read the actual code. The general structure is shown below. ```python def func name (...parameters...): ””” Usage:... Args:... Returns:... ””” ``` ```python def func name (...parameters...): ””” Usage:... Args:... Returns:... ””” ``` ... ... d) knobase.txt: This is different for every problem and is written by LLM Researcher. K. Generated Codes Comparison with other LHHs To further exemplify the algorithmic complexity differences mentioned above, we compare the MIS solver generated by BEAM and EoH. From Table V, we can find that the performance gap between BEAM and EoH widens significantly on harder instances (RB 800-1200), as shown in Table V. This divergence stems from EoH’s oversimplified crossover mechanism - it relies exclusively on uniform crossover (See Fig. 10). While this simplicity may leave more computational budget for RLSA local search on smaller instances (RBSmall), it fundamentally limits EoH’s ability to escape local optima on larger, more challenging instances (RB-Large). In contrast, BEAM’s more sophisticated evolutionary framework enables stronger capabilities, leading to consistently better performance as problem difficulty increases.

17

For TSP and CVRP, the conclusions are similar, with Fig. 11 illustrating the iteration curve of the generated algorithms (the curve is averaged over 5 generated algorithms). L. Generated Codes Introduction a) TSP - Traditional Benchmark: The best algorithm (Lisiting 1) features a two-stage computation process that first calculates edge utilities from normalized distances and then transforms them into penalties using configurable hyperparameters (ALPHA, BETA) and non-linear scaling. This algorithm is generated within a fixed budget, so it isn’t necessarily the best. b) BPP - Traditional Benchmark: The best algorithm (Listing 2) implements a three-phase, time-aware bin selection strategy for Bin Packing. It blends capacity-based and fit-based scoring, leverages a fast exact-fit shortcut, aggressively avoids overfilling bins via a lookahead penalty, and refines selections under tight time constraints. This algorithm is generated within a fixed budget, so it isn’t necessarily the best. Despite its good performance compared with its counterparts, we must note that BEAM tends to overcomplicate solutions for simple objectives - a tendency clearly reflected in code length. While EoH-generated solutions typically maintain concise implementations under 20 lines, BEAM’s output often exhibits unnecessary complexity. c) CAF - Traditional Benchmark: The best algorithm (Listing 3) dynamically adjusts exploration and exploitation priorities based on optimization progress, which is jointly characterized by budget consumption and solution quality. By blending standard and phase-aware Expected Improvement (EI) and scaling cost penalties according to phase, the method ensures robust and adaptive decision-making. The utility function further amplifies this adaptivity through exponentiation and diminishing sensitivity to cost over time. d) MIS with RLSA: The best algorithm (Listing 4) is overall a standard memetic algorithm, standing out by evolving entirely within the feasible independent set space—thanks to heuristic initialization, conflict-free crossover, and RLSAdriven local search—which accelerates convergence and ensures consistently valid, high-quality solutions. A comparison with EoH is given in Fig. 10. e) MIS with KaHIP & ARW: The best algorithm (Listing 5) utilizes KaHIP [66] in the initialization stage by isolating populations per partition and enabling occasional inter-block exchange, which is the best usage of KaHIP BEAM has found. Another interesting finding is that this algorithm, which outperforms KaMIS, relies on a simple uniform crossover—unlike KaMIS [13], which uses KaHIP specifically in the crossover stage. To further investigate, we manually replaced the uniform crossover with KaHIP-based crossover, mimicking KaMIS’s approach. Surprisingly, this change led to worse performance, suggesting that utilizing KaHIP in the crossover stage doesn’t necessarily boost the performance. f) CVRP with Split & LS: The best algorithm (Listing 6) combines 4 initialization strategies to create a diverse initial population of solutions. It features an adaptive perturbation mechanism [67] that employs 4 different mutation

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

18

Fig. 10: EoH vs. BEAM: MIS solver with RLSA.

strategies with problem-size-dependent intensity for effective exploration. The implementation also incorporates periodic intensification to refine good solutions. It’s worth noting that this codes serve as a good example of MCTS’s strength since from the comment in func_1 we can know that LLM only plans to realize two mutation strategies (swap and reverse) in exterior structure evolution. When realizing the function in the interior layer, LLM expands the strategy sets with two more complicated strategies (shift and scramble). Similarly, func_2 provides a correct evaluation function, which is also meaningful since we find that singlelayered EoH may even output the wrong evaluation function. g) TSP with EDM: The best algorithm (Listing 7) implements a well-designed ACO framework. Its core strength lies

in the balanced integration of pheromone-guided exploration and adaptive 2-opt refinement, enabling robust performance across diverse problem scales. The algorithm preserves solution quality via elite-preservation mechanisms and on-demand local optimization h) BBOB: The best algorithm (Listing 8) intelligently combines differential evolution, CMA-ES, and PSO in a staged optimization framework. Its key advantage lies in the dynamic allocation of computational budgets to each method based on their complementary strengths - DE for broad exploration, CMA-ES for precise local refinement, and PSO for final polishing. It automatically adjusts critical parameters during optimization, delivering robust performance across diverse continuous optimization landscapes. A smart restart mecha-

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

19

Fig. 11: Generated Algorithm Iteration Curve.

nism further enhances solution quality.

penalties = init_penalty * penalty_adjustment return penalties

M. Codes generated by BEAM Listing 1: BEAM-generated penalty function for TSP-GLS. ALPHA = 0.3 BETA = 1.5 INIT_PENALTY = 1 import numpy as np def heuristic(distance_matrix: np.ndarray) -> np.ndarray: # The ‘heuristic‘ function takes as input a distance ,→ matrix, and returns prior indicators of how bad it ,→ is to include each edge in a solution. The return is ,→ of the same shape as the input. # Step 1: Calculate initial edge utilities based on ,→ distance matrix edge_utilities = func_1(distance_matrix) # Step 2: Compute penalties based on edge utilities and ,→ current penalties penalties = func_2(edge_utilities, ALPHA, BETA, ,→ INIT_PENALTY) return penalties def func_1(distance_matrix: np.ndarray) -> np.ndarray: # Purpose: Calculate initial edge utilities based on ,→ distance matrix (higher distance = higher utility) # Normalize the distance matrix to [0,1] range for ,→ utility calculation max_dist = np.max(distance_matrix) if max_dist > 0: normalized_dist = distance_matrix / max_dist else: normalized_dist = distance_matrix return normalized_dist def func_2(edge_utilities: np.ndarray, alpha: float, beta: ,→ float, init_penalty: float) -> np.ndarray: # Purpose: Compute penalties for edges based on their ,→ utilities and hyperparameters # Calculate penalty adjustment factor utility_factor = alpha * edge_utilities # Apply non-linear transformation using beta exponent penalty_adjustment = np.power(utility_factor, beta) # Compute final penalties by scaling with initial ,→ penalty

Listing 2: BEAM-generated BPP priority function. MAX_TIME = 2 FITNESS_WEIGHT = 0.5514079756555896 CAPACITY_WEIGHT = 0.45840214751046593 STABILITY_WEIGHT = 0.10166942364627612 FILL_THRESHOLD = 0.7993217820159937 TIME_CHECK_INTERVAL = 0.0021170327307504515 LOOKAHEAD_PENALTY = 0.4432174448805559 CORE_PHASE_RATIO = 0.6278250358826285 import numpy as np import time def heuristic(item: float, bins_remain_cap: np.ndarray) -> ,→ np.ndarray: """ Hybrid heuristic combining the best elements from both ,→ approaches: 1. Maintains three-phase structure but with improved ,→ time allocation 2. Uses capacity-based scoring inspired by elitist code 3. Incorporates more aggressive fill threshold from ,→ elitist version 4. Optimized time checks and weight distribution """ start_time = time.time() # Phase 0: Fast exact fit check (immediate return if ,→ found) if time.time() - start_time > MAX_TIME: return np.zeros_like(bins_remain_cap) exact_fit_mask = (bins_remain_cap == item) if np.any(exact_fit_mask): return np.where(exact_fit_mask, np.inf, -np.inf) # Initialize scores with capacity validation valid_bins = bins_remain_cap >= item priority_scores = np.where(valid_bins, 0.0, -np.inf) # Phase 1: Core calculations (time-constrained) if time.time() - start_time < MAX_TIME * ,→ CORE_PHASE_RATIO: # Parallel score calculations with frequent time ,→ checks capacity_scores = func_1(bins_remain_cap, item) if time.time() - start_time > MAX_TIME:

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

return np.where(valid_bins, capacity_scores, ,→ np.inf) fitness_scores = func_2(bins_remain_cap, item) if time.time() - start_time > MAX_TIME: combined = FITNESS_WEIGHT * fitness_scores + ,→ CAPACITY_WEIGHT * capacity_scores return np.where(valid_bins, combined, -np.inf) # Combine scores with weights priority_scores = np.where( valid_bins, FITNESS_WEIGHT * fitness_scores + ,→ CAPACITY_WEIGHT * capacity_scores, -np.inf ) # Apply adaptive penalty (stronger than original) nearly_full = (bins_remain_cap / (bins_remain_cap + ,→ item)) > FILL_THRESHOLD priority_scores[nearly_full] *= LOOKAHEAD_PENALTY # Phase 2: Strategic refinement (if time permits) if time.time() - start_time < MAX_TIME * 0.9: priority_scores = func_3(priority_scores, ,→ bins_remain_cap, item)

20

best_x: torch.Tensor, best_y: int, test_x: torch.Tensor, mean_test_y: torch.Tensor, std_test_y: torch.Tensor, cost_test_y: torch.Tensor, budget_used: int, budget_total: int ) -> torch.Tensor: """ Optimized heuristic combining: 1. Dynamic phase calculation with budget and quality ,→ awareness 2. Balanced exploration-exploitation tradeoff 3. Phase-aware non-linear utility combination 4. Strict time constraints with frequent checks 5. Cost penalty that increases with phase """ start_time = time.time() # Calculate optimization phase considering both budget ,→ and solution quality phase = func_1(budget_used, budget_total, best_y, ,→ train_y, PHASE_SMOOTHING) if time.time() - start_time > MAX_TIME: return torch.zeros_like(mean_test_y)

return priority_scores def func_1(bins_remain_cap: np.ndarray, item: float) -> np. ,→ ndarray: # Purpose: Calculate capacity-based priority ( ,→ normalized remaining capacity) max_cap = np.max(bins_remain_cap) if max_cap == 0: return np.zeros_like(bins_remain_cap) return bins_remain_cap / max_cap def func_2(bins_remain_cap: np.ndarray, item: float) -> np. ,→ ndarray: # Purpose: Calculate fit-based priority (1 - abs( ,→ remaining_cap - item)/item) if item == 0: return np.zeros_like(bins_remain_cap) fit_quality = 1 - np.abs(bins_remain_cap - item) / item return np.maximum(0, fit_quality) # Ensure non,→ negative scores def func_3(priority_scores: np.ndarray, bins_remain_cap: np ,→ .ndarray, item: float) -> np.ndarray: # Purpose: Apply look-ahead adjustment considering ,→ potential future items # Penalize bins that would leave too little remaining ,→ capacity for typical future items remaining_after_packing = bins_remain_cap - item avg_item_size = item # Using current item as estimate future_fit_penalty = np.where( remaining_after_packing < avg_item_size * 0.5, 0.7, # Strong penalty if unlikely to fit another ,→ item 1.0 # No penalty if likely to fit another item ) return priority_scores * future_fit_penalty

Listing 3: BEAM-generated CAF. from heubase.caf import EI # This function is designed by ,→ BEAM from heubase.caf import phase_aware_EI # This function is ,→ designed by BEAM from heubase.caf import phase_aware_cost_scaling # This ,→ function is designed by BEAM MAX_TIME = 2 COST_DISCOUNT_FACTOR = 0.4506552556778862 IMPROVEMENT_BOOST = 3.2386305441223917 EXPLORATION_FACTOR = 0.3038707156753378 UTILITY_EXPONENT = 1.8344236311494773 PHASE_SMOOTHING = 0.6121269499903754 COST_PENALTY = 0.12516595027522082 import torch import time def heuristic(train_x: torch.Tensor, train_y: torch.Tensor,

# Get both base and boosted EI values base_ei = EI(mean_test_y, std_test_y, best_y) boosted_ei = phase_aware_EI( mean=mean_test_y, std=std_test_y, best_y=best_y, phase=phase, improvement_boost=IMPROVEMENT_BOOST ) if time.time() - start_time > MAX_TIME: return torch.zeros_like(mean_test_y) # Dynamic EI blending based on phase exploration_weight = EXPLORATION_FACTOR * (1 - phase) ei_values = (1 - exploration_weight) * boosted_ei + ,→ exploration_weight * base_ei # Get phase-sensitive cost scaling with additional ,→ penalty scaled_costs = phase_aware_cost_scaling( cost=cost_test_y, budget_used=budget_used, budget_total=budget_total, phase=phase, cost_discount_factor=COST_DISCOUNT_FACTOR ) * (1 + COST_PENALTY * phase) if time.time() - start_time > MAX_TIME: return torch.zeros_like(mean_test_y) # Phase-adaptive utility combination with exponent ,→ control utility = func_2(ei_values, scaled_costs, phase, ,→ UTILITY_EXPONENT) return utility def func_1(budget_used: int, budget_total: int, best_y: ,→ float, train_y: torch.Tensor, smoothing: float) -> float ,→ : """ Purpose: Calculate comprehensive optimization phase considering: 1. Budget consumption ratio (linear) 2. Solution quality improvement (non-linear) 3. Smooth transitions between phases Returns normalized phase value [0,1] """ # Budget-based phase component budget_phase = min(budget_used / budget_total, 1.0) # Quality-based phase component (normalized improvement ,→ ) min_y = train_y.min().item() max_y = train_y.max().item() if max_y != min_y: quality_phase = (best_y - min_y) / (max_y - min_y)

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

else: quality_phase = 0.0 # Combined phase with smoothing combined_phase = smoothing * budget_phase + (1 ,→ smoothing) * quality_phase return min(max(combined_phase, 0.0), 1.0) def func_2(ei_values: torch.Tensor, scaled_costs: torch. ,→ Tensor, phase: float, exponent: float) -> torch.Tensor: """ Purpose: Advanced utility function that: 1. Uses exponential scaling for non-linearity 2. Adapts cost sensitivity based on phase 3. Maintains numerical stability 4. Incorporates diminishing returns on high-cost ,→ evaluations """ # Phase-adaptive cost sensitivity cost_sensitivity = 1.0 - 0.5 * phase # Reduces cost ,→ sensitivity as optimization progresses # Numerically stable utility calculation with exponent eps = 1e-8 utility = (ei_values + eps).pow(exponent) / ( ,→ scaled_costs + eps).pow(cost_sensitivity) return utility

Listing 4: BEAM-generated MIS solver with RLSA. from rlsa import rlsa import torch from torch import Tensor import numpy as np import time MAX_TIME = 120 # Manullay set by us RLSA_ITERATIONS = 150 # Manullay set by us POPULATION_SIZE = 47 MUTATION_RATE = 0.08673770360907226 CROSSOVER_RATE = 0.8442911280296675 TOURNAMENT_SIZE = 5

def func_1(num_nodes: int, edge_index: np.ndarray) -> np. ,→ ndarray: # Purpose: Generate a random valid initial solution ( ,→ independent set) solution = np.zeros(num_nodes, dtype=int) adj_list = [[] for _ in range(num_nodes)] for i in range(edge_index.shape[1]): u, v = edge_index[0, i], edge_index[1, i] adj_list[u].append(v) adj_list[v].append(u) nodes = np.random.permutation(num_nodes) for node in nodes: if not any(solution[neighbor] for neighbor in ,→ adj_list[node]): solution[node] = 1 return solution def func_2(parent1: np.ndarray, parent2: np.ndarray, ,→ num_nodes: int, edge_index: np.ndarray) -> np.ndarray: # Purpose: Crossover two parents via union and conflict ,→ resolution parent1 = np.asarray(parent1, dtype=int) parent2 = np.asarray(parent2, dtype=int) child = np.zeros(num_nodes, dtype=int) adj_list = [[] for _ in range(num_nodes)] for i in range(edge_index.shape[1]): u, v = edge_index[0, i], edge_index[1, i] adj_list[u].append(v) adj_list[v].append(u) # Merge parents using a logical OR. candidates = np.where(parent1.astype(bool) | parent2. ,→ astype(bool))[0] np.random.shuffle(candidates) selected = [] for node in candidates: if not any(child[neighbor] for neighbor in adj_list ,→ [node]):

21

selected.append(node) child[node] = 1 return child def func_3(solution: np.ndarray, num_nodes: int, edge_index ,→ : np.ndarray, mutation_rate: float) -> np.ndarray: # Purpose: Mutate by flipping nodes while maintaining ,→ validity mutated = np.asarray(solution, dtype=int).copy() adj_list = [[] for _ in range(num_nodes)] for i in range(edge_index.shape[1]): u, v = edge_index[0, i], edge_index[1, i] adj_list[u].append(v) adj_list[v].append(u) for node in range(num_nodes): if np.random.rand() < mutation_rate: if mutated[node] == 1: mutated[node] = 0 else: if all(mutated[neighbor] == 0 for neighbor ,→ in adj_list[node]): mutated[node] = 1 return mutated def func_4(graph, x: Tensor, penalty_coeff: float) -> Tuple ,→ [Tensor, Tensor]: x_uq = x.unsqueeze(1) energy_term1 = torch.sum(x, dim=1) energy_term2 = torch.sum((torch.matmul(x_uq, graph) * ,→ x_uq).squeeze(1), 1) energy = -energy_term1 + penalty_coeff * energy_term2 grad_term1 = torch.ones_like(x) grad_term2 = penalty_coeff * torch.matmul(graph, x. ,→ unsqueeze(-1)).squeeze(-1) grad = -grad_term1 + grad_term2 return energy, grad def heuristic(num_nodes: int, edge_index: np.ndarray) -> np ,→ .ndarray: start_time = time.time() best_solution = np.zeros(num_nodes, dtype=int) best_fitness = 0 # Initialize population population = [func_1(num_nodes, edge_index) for _ in ,→ range(POPULATION_SIZE)] fitness = [ind.sum() for ind in population] best_idx = np.argmax(fitness) best_solution, best_fitness = population[best_idx].copy ,→ (), fitness[best_idx] print(f"Initial best: {best_fitness}") while time.time() - start_time < MAX_TIME: # Parent selection (tournament) parents = [] for _ in range(POPULATION_SIZE): candidates = np.random.choice(POPULATION_SIZE, ,→ TOURNAMENT_SIZE, replace=False) best = candidates[np.argmax([fitness[c] for c ,→ in candidates])] parents.append(population[best]) # Crossover and mutation offspring = [] for i in range(0, POPULATION_SIZE, 2): p1, p2 = parents[i], parents[min(i+1, ,→ POPULATION_SIZE-1)] if np.random.rand() < CROSSOVER_RATE: c1 = func_2(p1, p2, num_nodes, edge_index) c2 = func_2(p2, p1, num_nodes, edge_index) else: c1, c2 = p1.copy(), p2.copy() c1 = func_3(c1, num_nodes, edge_index, ,→ MUTATION_RATE) c2 = func_3(c2, num_nodes, edge_index, ,→ MUTATION_RATE) c1 = rlsa(num_nodes, edge_index, c1, ,→ RLSA_ITERATIONS, True, func_4) c2 = rlsa(num_nodes, edge_index, c2, ,→ RLSA_ITERATIONS, True, func_4) offspring.extend([c1, c2]) # Evaluate offspring offspring_fitness = [c.sum() for c in offspring]

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

combined_pop = population + offspring combined_fit = fitness + offspring_fitness # Survival selection sorted_idx = np.argsort(combined_fit)[::-1][: ,→ POPULATION_SIZE] population = [combined_pop[i] for i in sorted_idx] fitness = [combined_fit[i] for i in sorted_idx] # Update best current_best = np.max(fitness) if current_best > best_fitness: best_fitness = current_best best_solution = population[np.argmax(fitness)]. ,→ copy() print(f"Best after iteration: {current_best}") # Time check if time.time() - start_time >= MAX_TIME: break return best_solution

Listing 5: BEAM-generated MIS solver with KaHIP & ARW. import numpy as np import time import random from heubase.kahip import node_separator from heubase.arw import arw MAX_TIME = 60 POPULATION_SIZE = 5 BLOCK_SIZE = 4 IMBALANCE = 0.1 SEED = 42 MUTATION_RATIO = 0.05 EXCHANGE_RATE = 0.2 # Calibration doesn’t offer an increase in performance so ,→ we adopt the initial setting. def func_1(xadj: np.ndarray, adjncy: np.ndarray, solution: ,→ np.ndarray) -> np.ndarray: num_nodes = len(solution) for i in range(num_nodes): if solution[i] == 1: for j in range(xadj[i], xadj[i+1]): neighbor = adjncy[j] if solution[neighbor] == 1: solution[neighbor] = 0 return solution def func_2(parent1: np.ndarray, parent2: np.ndarray, xadj: ,→ np.ndarray, adjncy: np.ndarray) -> tuple[np.ndarray, np. ,→ ndarray]: # Purpose: Uniform crossover with mask generation mask = np.random.randint(0, 2, parent1.size).astype( ,→ bool) child1 = np.where(mask, parent1, parent2) child2 = np.where(mask, parent2, parent1) return child1, child2 def heuristic(num_nodes, edge_index, xadj, adjncy): start_time = time.time() best_solution = np.zeros(num_nodes, dtype=int) best_size = 0 # Initialize BLOCK_SIZE separate populations based on ,→ graph partitions _, parts = node_separator(xadj, adjncy, num_blocks= ,→ BLOCK_SIZE, imbalance=IMBALANCE, seed=SEED, mode=0) populations = [] for block in range(BLOCK_SIZE): population = [] for _ in range(POPULATION_SIZE): ind = np.zeros(num_nodes, dtype=int) ind[parts == block] = np.random.randint(0, 2, ,→ size=np.sum(parts == block)) population.append(func_1(xadj, adjncy, ind)) populations.append(population) while time.time() - start_time < 60: # Evolve each population separately for block in range(BLOCK_SIZE): # Evaluate

22

fitness = [np.sum(ind) for ind in populations[ ,→ block]] best_idx = np.argmax(fitness) if fitness[best_idx] > best_size: best_solution = populations[block][best_idx ,→ ].copy() best_size = fitness[best_idx] print(f"Best size: {best_size}") # Selection new_pop = [] for _ in range(POPULATION_SIZE): a, b = random.sample(range(POPULATION_SIZE) ,→ , 2) winner = a if fitness[a] > fitness[b] else ,→ b new_pop.append(populations[block][winner]. ,→ copy()) # Crossover for i in range(0, POPULATION_SIZE, 2): if i+1 >= POPULATION_SIZE: break child1, child2 = func_2(new_pop[i], new_pop ,→ [i+1], xadj, adjncy) new_pop[i] = func_1(xadj, adjncy, child1) new_pop[i+1] = func_1(xadj, adjncy, child2) # Mutation for ind in new_pop: for _ in range(int(num_nodes * ,→ MUTATION_RATIO)): idx = random.randint(0, num_nodes-1) ind[idx] = 1 - ind[idx] func_1(xadj, adjncy, ind) # Local search on best new_pop[0] = arw(xadj, adjncy, new_pop[0]) populations[block] = new_pop # Periodically exchange individuals between ,→ populations if random.random() < EXCHANGE_RATE: src, dest = random.sample(range(BLOCK_SIZE), 2) idx = random.randint(0, BLOCK_SIZE) populations[dest][idx] = populations[src][idx]. ,→ copy() return best_solution

Listing 6: BEAM-generated CVRP solver with Split & Local Search. from heubase.hgs import split from heubase.hgs import LS_Valid from heubase.hgs import LS_Invalid from heubase.cvrp import sweep_init # This function is ,→ designed by BEAM import numpy as np from typing import Tuple MAX_TIME = 300 INITIAL_POOL_SIZE = 6 PERTURB_STRENGTH = 0.15 LS_INTENSIFY_PROB = 0.5 import time import random def func_1(perm: np.ndarray) -> np.ndarray: # Purpose: Perform adaptive permutation perturbation ,→ using swap and reverse mutations n = len(perm) perturbed = perm.copy() # Determine perturbation strength based on problem size k = max(1, int(n * PERTURB_STRENGTH)) # Randomly choose between different perturbation ,→ strategies strategy = np.random.choice([’swap’, ’reverse’, ’shift’ ,→ , ’scramble’]) if strategy == ’swap’: # Perform k random swaps

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

for _ in range(k): i, j = np.random.choice(n, 2, replace=False) perturbed[i], perturbed[j] = perturbed[j], ,→ perturbed[i] elif strategy == ’reverse’: # Reverse a random subsequence i = np.random.randint(0, n - k + 1) perturbed[i:i+k] = perturbed[i:i+k][::-1] elif strategy == ’shift’: # Shift a random subsequence to a new position i = np.random.randint(0, n - k + 1) j = np.random.randint(0, n - k + 1) while abs(i - j) < k: # Ensure meaningful shift j = np.random.randint(0, n - k + 1) segment = perturbed[i:i+k] remaining = np.delete(perturbed, slice(i, i+k)) insert_pos = j if j < i else j - k perturbed = np.insert(remaining, insert_pos, ,→ segment) elif strategy == ’scramble’: # Scramble a random subsequence i = np.random.randint(0, n - k + 1) segment = perturbed[i:i+k] np.random.shuffle(segment) perturbed[i:i+k] = segment return perturbed def func_2(solution: np.ndarray, dist_matrix: np.ndarray) ,→ -> float: # Purpose: Calculate total route distance using ,→ precomputed distance matrix total_distance = 0.0 prev_node = 0 # Start at depot for node in solution[1:]: # Skip first depot (already ,→ accounted for in prev_node) total_distance += dist_matrix[prev_node, node] prev_node = node # Add return to depot from last node total_distance += dist_matrix[prev_node, 0] return total_distance def heuristic( nbClients: int, nbVehicles: int, capacity: float, depot_coord: np.ndarray, nodes_coord: np.ndarray, demands: np.ndarray, ) -> np.ndarray: start_time = time.time() best_sol = None best_cost = float(’inf’) # Precompute distance matrix all_nodes = np.vstack([depot_coord, nodes_coord]) dist_matrix = np.sqrt(((all_nodes[:, np.newaxis] ,→ all_nodes)**2).sum(axis=2)) # Generate diverse initial permutations initial_perms = [] # 1. Nearest neighbor heuristic if time.time() - start_time < MAX_TIME and nbClients > ,→ 0: current = 0 unvisited = set(range(1, nbClients+1)) nn_perm = [] while unvisited: nearest = min(unvisited, key=lambda x: ,→ dist_matrix[current, x]) nn_perm.append(nearest) unvisited.remove(nearest) current = nearest initial_perms.append(np.array(nn_perm)) # 2. Farthest insertion heuristic if time.time() - start_time < MAX_TIME and nbClients > ,→ 0: farthest = np.argmax(dist_matrix[0, 1:]) + 1

23

ff_perm = [farthest] unvisited = set(range(1, nbClients+1)) - {farthest} while unvisited and time.time() - start_time < ,→ MAX_TIME: candidate = max(unvisited, key=lambda x: min( ,→ dist_matrix[x][y] for y in ff_perm)) ff_perm.append(candidate) unvisited.remove(candidate) initial_perms.append(np.array(ff_perm)) # 3. Sweep algorithm if time.time() - start_time < MAX_TIME and nbClients > ,→ 0: initial_perms.append(sweep_init(nodes_coord)) # 4. Random permutations while len(initial_perms) < INITIAL_POOL_SIZE and time. ,→ time() - start_time < MAX_TIME: initial_perms.append(np.random.permutation( ,→ nbClients) + 1) # Evaluate initial solutions for perm in initial_perms: if time.time() - start_time >= MAX_TIME: break solution, valid = split(nbClients, nbVehicles, ,→ capacity, depot_coord, nodes_coord, demands, ,→ perm) if valid: improved_sol, improved_valid = LS_Valid( ,→ nbClients, nbVehicles, capacity, depot_coord ,→ , nodes_coord, demands, solution) current_cost = func_2(improved_sol, dist_matrix ,→ ) if improved_valid else float(’inf’) if current_cost < best_cost: best_sol = improved_sol best_cost = current_cost print(f"Initial best: {best_cost}") # Main optimization loop while time.time() - start_time < MAX_TIME: if best_sol is None: # Fallback initialization perm = np.random.permutation(nbClients) + 1 solution, valid = split(nbClients, nbVehicles, ,→ capacity, depot_coord, nodes_coord, demands, ,→ perm) if valid: best_sol, best_cost = solution, func_2( ,→ solution, dist_matrix) continue # Perturb current best permutation current_perm = best_sol[best_sol != 0].astype(int) perturbed_perm = func_1(current_perm) # Split and improve new_sol, new_valid = split(nbClients, nbVehicles, ,→ capacity, depot_coord, nodes_coord, demands, ,→ perturbed_perm) if new_valid: improved_sol, improved_valid = LS_Valid( ,→ nbClients, nbVehicles, capacity, depot_coord ,→ , nodes_coord, demands, new_sol) else: improved_sol, improved_valid = LS_Invalid( ,→ nbClients, nbVehicles, capacity, depot_coord ,→ , nodes_coord, demands, new_sol) # Evaluate and update if improved_valid: current_cost = func_2(improved_sol, dist_matrix ,→ ) if current_cost < best_cost: best_sol = improved_sol best_cost = current_cost print(f"Iteration best: {best_cost}") # Periodic intensification if np.random.rand() < LS_INTENSIFY_PROB and ,→ best_sol is not None: intensified_sol, intens_valid = LS_Valid( ,→ nbClients, nbVehicles, capacity, depot_coord ,→ , nodes_coord, demands, best_sol) if intens_valid: intens_cost = func_2(intensified_sol, ,→ dist_matrix)

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

if intens_cost < best_cost: best_sol = intensified_sol best_cost = intens_cost print(f"Intensification best: { ,→ best_cost}") # Final validity check if best_sol is None: perm = initial_perms[0] if initial_perms else np. ,→ arange(1, nbClients+1) best_sol, _ = split(nbClients, nbVehicles, capacity ,→ , depot_coord, nodes_coord, demands, perm) return best_sol.astype(int) if best_sol is not None ,→ else np.array([0], dtype=int)

Listing 7: BEAM-generated TSP solver with EDM import numpy as np from scipy.spatial import distance_matrix from heubase.tsp_aco import compute_tsp_edge_heuristics from heubase import two_opt_local_search # This function is ,→ designed by BEAM NUM_ANTS = 100 EVAPORATION_RATE = 0.1 ALPHA = 1.0 BETA = 2.0 INIT_PHEROMONE = 0.1 Q = 100.0 ELITE_FACTOR = 2.0 MAX_TIME = 30 import numpy as np import time def heuristic(node_coor: np.ndarray) -> np.ndarray: n = node_coor.shape[0] dist_matrix = np.sqrt(((node_coor[:, None] - node_coor) ,→ **2).sum(axis=2)) np.fill_diagonal(dist_matrix, 1.0) heuristic_matrix = compute_tsp_edge_heuristics( ,→ dist_matrix) pheromone = np.full((n, n), INIT_PHEROMONE) best_path = None best_length = np.inf start_time = time.time() while time.time() - start_time < MAX_TIME: paths = [] lengths = [] # Generate ant paths for _ in range(NUM_ANTS): current = np.random.randint(n) path = [current] unvisited = set(range(n)) - {current} while unvisited: current_node = path[-1] next_node = func_1( current_node, list(unvisited), pheromone[current_node], ,→ heuristic_matrix[current_node], ALPHA, BETA ) path.append(next_node) unvisited.remove(next_node) # Apply local optimization optimized_path = two_opt_local_search(np.array( ,→ path), dist_matrix) if len(np.unique(optimized_path)) == n: path = optimized_path # Calculate open TSP length length = dist_matrix[path[:-1], path[1:]].sum() paths.append(path) lengths.append(length) # Update best solution if length < best_length: best_path = np.array(path) best_length = length # Pheromone update pheromone = func_2(

24

pheromone, paths, lengths, best_path, best_length, EVAPORATION_RATE, Q, ELITE_FACTOR ) # print(f"Best length: {best_length}, Path: { ,→ best_path}") return best_path.astype(int) if best_path is not None ,→ else np.arange(n) def func_1(current_node: int, unvisited: list, ,→ pheromone_row: np.ndarray, heuristic_row: np.ndarray, ,→ alpha: float, beta: float) -> int: # Purpose: Select next node using probabilistic rule ( ,→ pheromoneˆalpha * heuristicˆbeta) unvisited = np.array(unvisited) pheromone = pheromone_row[unvisited] heuristic = heuristic_row[unvisited] probabilities = (pheromone ** alpha) * (heuristic ** ,→ beta) probabilities /= probabilities.sum() return np.random.choice(unvisited, p=probabilities) def func_2(pheromone: np.ndarray, paths: list, lengths: ,→ list, best_path: np.ndarray, best_length: float, ,→ evaporation_rate: float, q: float, elite_factor: float) ,→ -> np.ndarray: # Purpose: Update pheromone with evaporation, ant ,→ deposits, and elite reinforcement # Evaporation pheromone *= (1 - evaporation_rate) # Ant deposits for path, length in zip(paths, lengths): delta = q / length for i in range(len(path)-1): u, v = path[i], path[i+1] pheromone[u, v] += delta pheromone[v, u] += delta # Elite reinforcement if best_path is not None: delta_elite = elite_factor * q / best_length for i in range(len(best_path)-1): u, v = best_path[i], best_path[i+1] pheromone[u, v] += delta_elite pheromone[v, u] += delta_elite return pheromone

Listing 8: BEAM-generated BBOB. import numpy as np from typing import Tuple, Callable import cma # pip-installed with the help of LLM Researcher BOUND = 5.12 P_DE = 0.6 P_CMA = 0.3 P_PSO = 0.1 POP_SIZE = 50 F_MIN = 0.4 F_MAX = 0.9 CR0 = 0.9 PSO_N = 10 W0 = 0.9 W1 = 0.4 C1 = 2.0 C2 = 2.0

def heuristic(problem_func: Callable, dimension: int, fopt: ,→ float, budget: int = 20000): bounds = (-BOUND, BOUND) evals = 0 bud_de = int(budget * P_DE) bud_cma = int(budget * P_CMA) bud_pso = budget - bud_de - bud_cma pop = bounds[0] + (bounds[1]-bounds[0]) * func_2(( ,→ POP_SIZE, dimension)) fit = np.array([problem_func(ind) for ind in pop]); ,→ evals += POP_SIZE

IEEE TRANSACTIONS ON NEURAL NETWORKS AND LEARNING SYSTEMS, VOL. XX, NO. X, APRIL 2026

best_idx = np.argmin(fit) best_x, best_f = pop[best_idx].copy(), fit[best_idx] iters_de = bud_de // POP_SIZE for t in range(1, iters_de+1): F = F_MIN + (F_MAX - F_MIN)*(1 - t/iters_de) CR = CR0 * np.exp(-3*t/iters_de) for i in range(POP_SIZE): idxs = [j for j in range(POP_SIZE) if j!=i] a,b,c = pop[np.random.choice(idxs,3,replace= ,→ False)] # current-to-best/1 mutant = pop[i] + F*(best_x-pop[i]) + F*(a-b) mutant = np.clip(mutant, bounds[0], bounds[1]) mask = np.random.rand(dimension) < CR if not mask.any(): mask[np.random.randint( ,→ dimension)] = True trial = np.where(mask, mutant, pop[i]) fv = problem_func(trial); evals+=1 if fv < fit[i]: pop[i], fit[i] = trial, fv if fv < best_f: best_x, best_f = trial.copy(), fv if t % 50 == 0: x2, f2 = func_1(best_x, problem_func, bounds) if f2 < best_f: best_x, best_f = x2.copy(), f2 sigma0 = np.std(pop, axis=0).mean() + 1e-8 es = cma.CMAEvolutionStrategy(best_x.tolist(), sigma0, ,→ {’popsize’:20}) while not es.stop() and evals < bud_de+bud_cma: X = es.ask() Fs = [problem_func(x) for x in X]; evals += len(X) es.tell(X, Fs) sol = np.array(es.result.xbest) fsol = problem_func(sol); evals+=1 if fsol < best_f: best_x, best_f = sol.copy(), fsol pos = best_x + 0.1 * np.random.randn(PSO_N, dimension) vel = np.zeros_like(pos) pbest = pos.copy() pfit = np.array([problem_func(x) for x in pos]); evals ,→ +=PSO_N gbest, gfit = best_x.copy(), best_f for k in range(bud_pso//PSO_N): w = W0 + (W1-W0)*(k/(bud_pso//PSO_N)) for i in range(PSO_N): r1, r2 = np.random.rand(dimension), np.random. ,→ rand(dimension) vel[i] = w*vel[i] + C1*r1*(pbest[i]-pos[i]) + ,→ C2*r2*(gbest-pos[i]) pos[i] = np.clip(pos[i] + vel[i], bounds[0], ,→ bounds[1]) fv = problem_func(pos[i]); evals+=1 if fv < pfit[i]: pbest[i], pfit[i] = pos[i].copy(), fv if fv < gfit: gbest, gfit = pos[i].copy(), fv best_x, best_f = gbest, gfit x3, f3 = func_1(best_x, problem_func, bounds, trials ,→ =20) if f3 < best_f: best_x, best_f = x3, f3 gap = best_f - fopt return best_x.tolist(), round(best_f,6), round(gap,6)

def func_1(x: np.array, prob: Callable, bounds: float, ,→ trials: int = 10) -> float: n = len(x) H = np.random.choice([1, -1], size=(trials, n)) cand = x + 1e-2 * H cand = np.clip(cand, bounds[0], bounds[1]) vals = [prob(c) for c in cand] idx = np.argmin(vals) return cand[idx], vals[idx] def func_2(shape: Tuple, mu: float = 0.3, iter: int = 5) -> ,→ Tuple: z = np.random.rand(*shape) for _ in range(iter): z = np.where(z < mu, z/mu, (1 - z)/(1 - mu))

return z

25