Paper deep dive
Distribution Matching Evolutionary Algorithms for Rare Event Sampling
Yonatan Gideoni, Yarin Gal
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 10/6/2026, 10:21:20 AM
Summary
This paper introduces Distribution Matching Evolutionary Algorithms (DME), a class of search methods that sample from a global target distribution without updating model weights. By interpreting evolutionary algorithms as approximate Markov Chain Monte Carlo (MCMC), the authors derive DME, which uses a sampling-resampling process to achieve higher sample efficiency than existing methods for rare event sampling and discovery problems.
Entities (7)
Relation Signals (5)
Yonatan Gideoni → affiliatedwith → University of Oxford
confidence 99% · Affiliation: University of Oxford
Yarin Gal → affiliatedwith → University of Oxford
confidence 99% · Affiliation: University of Oxford
Distribution Matching Evolutionary Algorithms → samplesfrom → global target distribution
confidence 95% · DME, a class of search methods which sample from a global target distribution without updating weights.
Distribution Matching Evolutionary Algorithms → interpretsas → Markov Chain Monte Carlo
confidence 90% · In this paper, we interpret various evolutionary algorithms as approximate Markov Chain Monte Carlo... This interpretation allows developing Distribution Matching Evolutionary Algorithms (DME)
Distribution Matching Evolutionary Algorithms → outperforms → existing methods
confidence 85% · Empirically, DME has a higher sample efficiency than existing methods on problems requiring many samples to find a solution.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:A novel discovery is one which is both useful and surprising: a generative model's output is a useful discovery if it has a low probability of being generated (it's surprising) and a high reward (it's useful). Global optimization can directly increase the probability of sampling high rewards but typically requires updating model weights. Such gradient based optimization is expensive and bars using capable closed-source models. Instead, modern search methods for discovery sacrifice the global target, and use evolutionary algorithms with local reward maximizing objectives, permitting the search to focus only on high probability samples. In this paper, we interpret various evolutionary algorithms as approximate Markov Chain Monte Carlo, an optimization-free method to sample from complex distributions. This interpretation allows developing Distribution Matching Evolutionary Algorithms (DME), a class of search methods which sample from a global target distribution without updating weights. Empirically, DME has a higher sample efficiency than existing methods on problems requiring many samples to find a solution.
Tags
Links
- Source: https://arxiv.org/abs/2610.03833v1
- Canonical: https://arxiv.org/abs/2610.03833v1
Trouble viewing inline? Open PDF directly →
Full Text
82,586 characters extracted from source content.
Expand or collapse full text
Distribution Matching Evolutionary Algorithms for Rare Event Sampling Yonatan Gideoni Yarin Gal Affiliation: University of Oxford Email: yg@robots.ox.ac.uk Abstract A novel discovery is one which is both useful and surprising: a generative model’s output is a useful discovery if it has a low probability of being generated (it’s surprising) and a high reward (it’s useful). Global optimization can directly increase the probability of sampling high rewards but typically requires updating model weights. Such gradient based optimization is expensive and bars using capable closed-source models. Instead, modern search methods for discovery sacrifice the global target, and use evolutionary algorithms with local reward maximizing objectives, permitting the search to focus only on high probability samples. In this paper, we interpret various evolutionary algorithms as approximate Markov Chain Monte Carlo, an optimization-free method to sample from complex distributions. This interpretation allows developing Distribution Matching Evolutionary Algorithms (DME), a class of search methods which sample from a global target distribution without updating weights. Empirically, DME has a higher sample efficiency than existing methods on problems requiring many samples to find a solution.11 1 For code see https://github.com/YonatanGideoni/dme. 1 Introduction This paper is concerned with the problem of how to find novel discoveries using generative models. Specifically, what does it mean for a discovery to be novel? Intuitively, a novel discovery is one which has never been seen before and is in some sense useful. Whether a discovery has been seen can be judged based on prior work, while its utility can be measured using some verifier or reward.22 2 There are very interesting, harder cases where how to measure utility is unclear – although highly relevant, we do not discuss those here. When searching for a discovery using a generative model, if it is novel then the model will likely assign it a low generation probability, due to not seeing the discovery in its training data. A common way to leverage generative models to find discoveries is using evolutionary algorithms (Novikov et al., 2025), see Algorithm 1 for an example. Evolutionary algorithms search for high reward samples by mutating existing candidates and selecting a subset based on their reward. Generative models are useful for mutating samples as they tend to yield semantically meaningful mutations, e.g. changing a program in a way so it still compiles. A common selection operator is to keep only the highest-reward candidates, which maximizes the local reward. Algorithm 1 Example Evolutionary Algorithm Input: An initial parent x0x_0, a mutation operator q(x′|x)q(x |x), a selection operator S, a fitness function F(x)F(x), population size N, number of generations n for i=0…ni=0...n do Generate a new population of offspring x1′,x2′,…,xN′∼q(⋅|xi)x _1,x _2,...,x _N q(·|x_i) Evaluate offspring fitnesses F(x1′),F(x2′),…,F(xN′)F(x _1),F(x _2),...,F(x _N) Select one of the offsprings as the next parent xi+1x_i+1 using S Return the fittest offspring Local reward maximization methods rarely sample low probability outputs due to getting stuck in local modes of the distribution. A better alternative to local reward maximization would be searching for not a high reward solution, but for all high reward solutions (Bengio et al., 2021). Concretely, instead of directly maximizing the reward, a more suitable goal would be to sample from a target distribution that focuses its probability mass on high reward regions. A mode collapsed search would sample solely from a single high reward region and thus be penalized. We refer to algorithms that sample from a target distribution as distribution matching algorithms. Generative models are typically trained for distribution matching using variational inference, where their weights are updated to minimize a KL divergence between the generated and the target distributions (Kingma and Welling, 2013; Ho et al., 2020; Malkin et al., 2022b; Eijkelboom et al., 2024). For discovery problems, training a generative model has two main downsides. First, it bars using capable fixed weights closed-source models, and second, training is much more expensive and time consuming than sampling. Due to these two reasons, most discovery methods are gradient-free evolutionary algorithms. A wide class of alternative, optimization-free methods for sampling from a target distribution are based on Markov Chain Monte Carlo (MCMC). Although widely used in statistics (Brooks et al., 2011), MCMC is less prevalent in modern generative modelling due to requiring many steps to generate a single sample (Goodfellow et al., 2014; Bengio et al., 2021). For discovery however, methods typically anyway require many steps and once the discovery is made there is no need to generate additional samples. Building on this insight, in this article we develop a class of inference only distribution matching algorithms. First, we demonstrate how some evolutionary algorithms can be interpreted as implementing approximate MCMC. We then use this interpretation to derive Distribution Matching Evolutionary Algorithms (DMEs), inference only search methods which sample from a target distribution. Empirically, DMEs are found to be more sample efficient than many existing methods on discovery problems which require many samples to find a solution. 2 Preliminaries This section introduces and builds intuition for some relevant methods. First, we explain how sampling can be cast as optimization, and how sampling from a given target distribution allows performing global optimization. Then we discuss a relaxation of this optimal target, yielding the target distribution of Bengio et al. (2021). Afterwards, we present MCMC, the statistics method used here to sample from the relaxed target distribution. 2.1 Sampling as Optimization Sampling from a target distribution can be used to optimize a target function, with the first and arguably most well known method for doing so being simulated annealing (Kirkpatrick et al., 1983). In simulated annealing, one samples from πT(x)∝exp(−βf(x)) _T(x) (-β f(x) ), where f(x)f(x) is the function to minimize and β=1/Tβ=1/T is an inverse temperature. As T is decreased (“annealed”) to 00, or equivalently β→∞β→∞, sampling from πT _T yields xxs that minimize f(x)f(x). Simulated annealing is useful but limited, as often sampling from πT _T is in itself nontrivial. In our setting, in place of some f(x)f(x) we presume access to some bounded positive reward function, R(x)R(x). To find all high reward solutions, one would ideally like to sample x uniformly from the max-reward set S≔x|R(x)=maxx′R(x′)S \x|R(x)= _x R(x )\, which would mean sampling from the idealized distribution πS(x)≔x∈S|S| _S(x) 1\x∈ S\|S|, where x∈S1\x∈ S\ is 11 when x∈Sx∈ S and 00 otherwise. Sampling from πS _S would result in sampling from all max-reward solutions in the space and never sampling a suboptimal solution. In this sense πS _S is the optimal target distribution, although it is unattainable in practice. 2.2 Relaxing the Optimal Target There are two main difficulties with using πS _S as a target distribution when searching for discoveries. First, a priori the maximum reward is often unknown, e.g. if searching for new mathematical bounds the tightest possible bound is unknown in advance. Second, πS(x) _S(x) encodes no notion of similarity or progress towards a higher reward the search algorithm may use, as it gives no probability mass to all suboptimal samples (see Figure 1 (b)). A straightforward relaxation of πS _S is π(x)=R(x)β/Zπ(x)=R(x)^β/Z, where Z is some unknown normalizing constant. As β→∞β→∞ one has that π(x)→πS(x)π(x)→ _S(x), as all suboptimal samples will have vanishingly small probability mass relative to the max-reward set. A higher β will thus better approximate the ideal target but yields a harder to sample from distribution, while lower β yield easier target distributions at the cost of more mass being placed on suboptimal solutions (see Figure 1 (c)). Figure 1: Given some reward function (a) sampling uniformly over all high reward solutions requires sampling from πS _S (b). πS _S is difficult to sample from due to slightly suboptimal solutions having no probability, so it is easier to sample from a relaxation (c). If relaxed too far, the target distribution will be almost uniform, leading to many wasted samples (d). Historically, using π(x)∝R(x)βπ(x) R(x)^β as a target distribution was first introduced by Bengio et al. (2021) for GFlowNets, although variations exist in earlier works (Kirkpatrick et al., 1983). GFlowNets aim to solve the same problem – sampling diverse high reward solutions – using methods from reinforcement learning and variational inference (Malkin et al., 2022b). Unlike simulated annealing, for GFlowNets β is typically kept fixed. Both methods approximate the target distribution by updating a proposal distribution’s parameters; next, we introduce a method to generate samples without needing to update the proposal distribution. 2.3 Markov Chain Monte Carlo To introduce MCMC it is helpful to first understand a simpler sampling algorithm, rejection sampling. Rejection sampling is a method to sample from some target distribution π(x)π(x) given a proposal distribution q(x)q(x). q is needed as typically π can be evaluated at a given point but cannot be sampled directly. Rejection sampling samples from π by first sampling x∼q(⋅)x q(·) and accepting the sample with probability min(1,π(x)Mq(x)) (1, π(x)Mq(x) ), where the π(x)q(x) π(x)q(x) term is an importance sampling correction and M is a constant which is sufficiently large to ensure that Mq(x)≥π(x)Mq(x)≥π(x). The sample is necessarily accepted if π(x)≥Mq(x)π(x)≥ Mq(x) and is only sometimes accepted if π(x)<Mq(x)π(x)<Mq(x). Although simple, rejection sampling is very sample inefficient, as in high dimensional spaces there may be many regions where the target’s density is much less than the proposal’s. In these cases π(x)≪Mq(x)π(x) Mq(x), making the acceptance probability very small. MCMC can be thought of as a sequential, local form of rejection sampling. MCMC is more robust to the choice of the proposal distribution than rejection sampling due to taking local steps, at the cost of needing multiple steps to generate a single independent sample from the target. The most general form of MCMC is Metropolis-Hastings MCMC (Metropolis et al., 1953; Hastings, 1970), where instead of accepting or rejecting a proposed point a change to an existing point is accepted or rejected.33 3 Unless mentioned otherwise, when discussing MCMC we specifically refer to Metropolis-Hastings MCMC. In MCMC the proposal distribution is replaced by a transition kernel, q(x′|x)q(x |x), and a proposed x′x is accepted with probability min(1,π(x′)q(x|x′)π(x)q(x′|x)) (1, π(x )q(x|x )π(x)q(x |x) ). The q(x′|x)q(x |x) terms arise as a correction factor, so the sampling will be unbiased even if it is likelier to transition from x to x′x than from x′x to x. If the transition kernel is symmetric, where q(x′|x)=q(x|x′)q(x |x)=q(x|x ), the acceptance probability simplifies to min(1,π(x′)π(x)) (1, π(x )π(x) ). Note that now the transition kernel drops out, so a proposal x′x is always accepted if it has a higher target probability than the previous point x, and is only sometimes accepted if its target probability is lower. Algorithm 2 shows the full MCMC loop. Algorithm 2 Metropolis-Hastings MCMC Input: A starting point x0x_0, a (potentially unnormalized) target distribution π(x)π(x), a transition kernel q(x′|x)q(x |x), number of steps n for i=0…ni=0...n do Sample x′∼q(⋅|xi)x q(·|x_i) Calculate acceptance probability a=min(1,π(x′)q(xi|x′)π(xi)q(x′|xi))a= (1, π(x )q(x_i|x )π(x_i)q(x |x_i) ) With probability a, set xi+1←x′x_i+1← x , otherwise xi+1←xix_i+1← x_i Return all xix_i; if n is sufficiently large, they form a (correlated) set of samples from π 3 Evolutionary Algorithms as Approximate Markov Chain Monte Carlo At a high level, evolutionary algorithms generate samples and select which to keep for the next generation. MCMC’s pseudocode is thus similar to that of certain evolutionary algorithms. In this section we show how some evolutionary algorithms can be interpreted as implementing approximate versions of MCMC, with a given evolutionary algorithm implicitly defining a transition kernel and a target distribution. While this interpretation will help us derive a class of evolutionary algorithms well suited for discovery, to build intuition we first demonstrate it using arguably the two simplest evolutionary algorithms, random search and hill climbing evolutionary strategies. The following derivations clarify which approximations the MCMC interpretation requires and how other, more general evolutionary algorithms can be interpreted as sampling from various target distributions as well. In all cases, we assume the transition kernel uses a generative model p(x)p(x) which can be conditioned using a prior sample, giving p(x′|x)p(x |x). For example, this could be a large language model (LLM) which is asked to generate code and is conditioned on a previously generated program. For ease of notation, we omit the generative model’s (fixed) parameters θ, the sampling parameters (e.g. the temperature), and other constant factors from the conditioning. 3.1 Random Search In random search one samples IID from the generative model, conditioned on no prior samples. Thus, the transition kernel is q(x′|x)=p(x′)q(x |x)=p(x ). As random sampling always accepts the generated samples, its acceptance probability is always 11, requiring 1=π(x′)q(x|x′)π(x)q(x′|x)=π(x′)p(x)π(x)p(x′)1= π(x )q(x|x )π(x)q(x |x)= π(x )p(x)π(x)p(x ). This is satisfied if π(x)=p(x)π(x)=p(x), which makes intuitive sense – given no signal, random sampling’s target is simply to sample from the proposal distribution. Interpreting random search as MCMC thus requires no approximations. 3.2 Hill Climbing Evolutionary Strategies There are two main kinds of hill climbing evolutionary strategies, (1+N)(1+N) and (1,N)(1,N) strategies (Beyer and Schwefel, 2002). The first, (1+N)(1+N), has N offspring generated from the parent and picks the best out of the N offspring+the one parent as the parent for the next generation. The second strategy, (1,N)(1,N), picks the best out of the N offspring, regardless of the parent’s fitness. Although MCMC by default generates only a single sample per step, it is possible extending it to use N samples per step. We extend MCMC to several samples using a trick used in some particle methods (Andrieu et al., 2010), known as sampled-importance-resampling (SIR) (Rubin, 1987). First, conditioned on the previous sample x, hill climbing samples N candidates from the generative model, yielding xi′∼p(⋅|x)x _i p(·|x) for 1≤i≤N1≤ i≤ N. To simulate the best-candidate choice used in hill climbing, we soft sample the max-reward by resampling a x′x from the xi′i=1N\x _i\_i=1^N by picking xi′x _i with probability proportional to R(xi′)βR(x_i )^β. Given a sufficiently high beta, this simulates the hard-max in hill climbing. Altogether, the sequential sampling-resampling process effectively yields the transition kernel q(x′|x)∝p(x′|x)R(x′)βq(x |x) p(x |x)R(x )^β, with hill climbing requiring β→∞β→∞. This transition kernel is equally valid for both (1+N)(1+N) and (1,N)(1,N) strategies. However, each strategy has a different acceptance criteria – (1+N)(1+N) accepts a proposal x′x only if it is better than x, while (1,N)(1,N) accepts it regardless. Starting from the (1,N)(1,N) case, to always accept a sample requires the target distribution to fulfill the property 1=π(x′)q(x|x′)π(x)q(x′|x)=π(x′)p(x|x′)R(x)βπ(x)p(x′|x)R(x′)β1= π(x )q(x|x )π(x)q(x |x)= π(x )p(x|x )R(x)^βπ(x)p(x |x)R(x )^β. As for generic p(x′|x)p(x |x) the always-accept condition cannot be fulfilled, to satisfy it requires either approximating or assuming the generative model’s conditional distribution is symmetric, so p(x′|x)=p(x|x′)p(x |x)=p(x|x ). Given a symmetric generative conditional distribution, to always-accept a sample requires π(x)∝R(x)βπ(x) R(x)^β. Note that this is the distribution-matching target distribution proposed in section 2, albeit at the infinite-beta limit. For (1+N)(1+N), to reject a sample if R(x)>R(x′)R(x)>R(x ) requires π(x)R(x)β≫π(x′)R(x′)β π(x)R(x)^β π(x )R(x )^β. In the infinite-beta limit, this is fulfilled by setting π(x)∝R(x)kβπ(x) R(x)^kβ for some k>1k>1. Thus, to summarize, hill climbing evolutionary strategies can be interpreted as MCMC with: 1. A transition kernel q(x′|x)∝p(x′|x)R(x′)βq(x |x) p(x |x)R(x )^β, where β→∞β→∞ 2. Assuming the generative model’s conditional distribution is symmetric: p(x′|x)=p(x|x′)p(x |x)=p(x|x ) 3. A target distribution π(x)∝R(x)kβπ(x) R(x)^kβ, where k>1k>1 for (1+N)(1+N) strategies and k=1k=1 for (1,N)(1,N) 4 Distribution Matching Evolutionary Algorithms In the previous sections we discuss what is a suitable distribution matching objective for discovery and introduce MCMC, which is necessary for sampling from the target distribution using inference-only methods. Additionally, we show how some evolutionary algorithms can be interpreted as implementing approximate MCMC. In this section we bring these tools together to design a practical algorithm to sample from the target using a generative model. The following can be interpreted as deriving a form of hill climbing which a) does not require a symmetric transition kernel and b) samples from the proxy target distribution π(x)∝R(x)βπ(x) R(x)^β given finite beta. Given the target distribution π(x)∝R(x)βπ(x) R(x)^β, what should be the transition kernel? The optimal kernel is the target, where q(x′|x)=π(x′)q(x |x)=π(x ), but if it was possible to sample from the target then there would be no need to use MCMC in the first place. However, it is possible approximately sampling from the optimal transition kernel by sampling in batches using sampled-importance-resampling. To illustrate, given x1,x2,…,xN∼Ux_1,x_2,...,x_N U, where U denotes the uniform distribution over the space, one can use batched importance sampling to sample from π instead. By sampling x′x from x1,…,xN\x_1,...,x_N\ with probability ∝R(xi)β R(x_i)^β, one effectively generates a sample from π. The sampling-resampling process will better approximate sampling from π a) for larger N and b) the closer π is to U, with the notion of closeness being made precise later. Sampling uniformly over the space would typically lead to very low quality samples, as most x would have a very low reward. A much better initial proposal would be to use the generative model, so x1,…,xN∼p(⋅|x)x_1,...,x_N p(·|x). Due to sampling from p(⋅|x)p(·|x) and not U, x′x requires an importance sampling correction, so samples should be selected from x1,…,xN\x_1,...,x_N\ with probability ∝R(xi)βp(xi|x) R(x_i)^βp(x_i|x). This process yields the transition kernel q(x′|x)≈p(x′|x)R(x′)βp(x′|x)Z=R(x′)βZ=π(x′),q(x |x)≈ p(x |x) R(x )^βp(x |x)Z= R(x )^βZ=π(x ), (1) where the first p(x′|x)p(x |x) stems from sampling xi∼p(xi|x)x_i p(x_i|x) while the R(x′)βp(x′|x)Z R(x )^βp(x |x)Z term comes from sampling x′x from x1,…,xN\x_1,...,x_N\ with probability ∝R(xi)βp(xi|x) R(x_i)^βp(x_i|x). Note that as q(x′|x)≈π(x′)q(x |x)≈π(x ) the MCMC acceptance probability is approximately one. In practice, we assume the approximation is exact and drop the accept/reject step. Altogether, this transition kernel, target distribution, and sampling process result in Algorithm 3. Note that as we are chiefly interested in discovery problems, only the sample with the highest reward is returned. Algorithm 3 Distribution Matching Evolution Input: A starting point x0x_0, a (potentially unnormalized) target distribution π(x)π(x), a generative model p(x′|x)p(x |x), batch size N, number of steps n for i=0…ni=0...n do Sample x1′,…,xN′∼p(⋅|xi)x _1,...,x _N p(·|x_i) Define the importance weights wi≔R(xi′)βp(xi′|xi)w_i R(x_i )^βp(x _i|x_i) and self-normalized probabilities pi≔wi∑j=1Nwjp_i w_i _j=1^Nw_j Sample x′∼x1′,…,xN′x \x _1,...,x _N\, where xi′x_i is sampled with probability pip_i Set xi+1←x′x_i+1← x Return the generated sample with the highest reward 4.1 Quality of Approximation The sampling-resampling process DME uses only approximates sampling from the true target distribution due to using a finite batch size. To understand how well sampling from the target is approximated we derive the first order correction to q(x′|x)q(x |x) in the sampling batch size, 1N 1N. Broadly, it is clear that if the generative model exactly samples from the target, so p(x′|x)=π(x′)p(x |x)=π(x ), there is no approximation, and in some sense the closer these two distributions are the better is the approximation. The size of this error is given by the first order correction, with proofs and mild necessary assumptions being delegated to Appendix A. Theorem 1 (First order correction to q(x′|x)q(x |x), informal). Define the transition kernel q(x′|x)q(x |x) using the sampling-resampling process given in Algorithm 3. Then, under some mild conditions, one has that q(x′|x)=π(x′)(1+D(x′|x)N)+O(N−3/2),q(x |x)=π(x ) (1+ D(x |x)N )+O (N^-3/2 ), (2) where D(x′|x)≔∫π(y)2p(y|x)y−π(x′)p(x′|x)D(x |x) π(y)^2p(y|x)dy- π(x )p(x |x). Note that D(x′|x)=0D(x |x)=0 if π(x′)=p(x′|x)π(x )=p(x |x), as required for the sampling process to be unbiased. To get a measure for the correction D(x′|x)D(x |x)’s magnitude, one can measure the average error when sampling proposals x′x from p(⋅|x)p(·|x) instead of from π. Doing so yields the following corollary, which quantifies this error. Corollary 2 (First order correction as a χ2χ^2 divergence). One has that x′∼p(⋅|x)[D(x′|x)]=∫(π(x′)−p(x′|x)p(x′|x))2p(x′|x)dx′=χ2(π||p(⋅|x))E_x p(·|x)[D(x |x)]= ( π(x )-p(x |x)p(x |x) )^2p(x |x)dx =χ^2(π||p(·|x)) (3) Where χ2χ^2 is the chi-squared divergence. Specifically, note that χ2χ^2 has the form of a weighted mean-square error between the two distributions. Thus, there are three ways to improve the approximate sampling from the target distribution: 1. Use an easier target which is more similar to the generative model’s distribution. Typically, this would mean lowering β. 2. Use a larger batch size N when sampling proposals. 3. Use a better generative model, or improve the given model’s conditioning. For example, when the generative model is a language model, a better prompt would focus more probability mass on high-reward regions. MCMC’s sequential sampling compensates for the transition kernel being only approximately optimal, where a better approximation requires fewer steps to effectively sample the target. 4.2 Novelty, Heuristic Interpretations, and Proposal-Agnosticism To better understand DME it helps interpreting its behavior in some limits, in which cases it reduces to various heuristics. For example, in the β→∞β→∞ limit DME would always select the highest reward sample, akin to hill climbing. Where it differs from typical hill climbing is in how ties are resolved given several samples with the same reward. If R(x1)=R(x2)R(x_1)=R(x_2) but p(x1|x)<p(x2|x)p(x_1|x)<p(x_2|x), DME will with higher probability pick x1x_1 due to the inverse dependence on the probability in the importance weights. In this sense DME prefers novel solutions, where lower probability samples are always picked all else being equal. If additionally one of the high-reward samples has a much lower probability than all the others, then DME is equivalent to hill climbing with a pick-lowest-probability tie-breaker. More generally, for finite β, DME can be interpreted as a form of stochastic hill climbing. An evolutionary algorithm typically cannot be described independently of its mutation operator, as different kinds of mutations can make the same algorithm give very different results. For example, hill climbing with one model proposing candidates may pick approach A 30%30\% of the time and approach B 70%70\%, while for a different model the split may instead be 90%/10%90\%/10\%. Unlike typical evolutionary algorithms DME is approximately proposal-agnostic, as the importance sampling weights correct for the probability the generative model gives to different samples. Due to the approximation’s quality, in practice DME may still be much more sample efficient with some models than with others. 5 DME for Language Models While DME as described in the previous section is suitable for generic generative models, in this section we describe two tweaks that make it better suited for language models. The general algorithm remains the same, with the two tweaks resulting only in modifying the importance sampling weights. Other design choices, for language models or other generative models, would accordingly give a different valid distribution matching evolutionary algorithm. 5.1 Path Marginalization While there are several ways to arrive at the same answer, typically the answer and not the path that led to it is what matters. For example, different reasoning trajectories may result in a language model writing the same code. If x is the previous code and x′x is the new code, although only the probability of generating the final code is of interest, the model on the way also generates some reasoning trajectory, z, yielding p(x′,z|x)p(x ,z|x) and not p(x′|x)p(x |x). More generally, z can be any intermediate steps towards the final answer, e.g. a set of tool calls and their returned outputs. The probability of arriving at the sampled code given the prior code, regardless of what the intermediate steps are, can be approximated as p(x′|x)=∫p(x′|z,x)p(z|x)dz⏟z∼p(⋅|x)[p(x′|z,x)]≈p(x′|z,x),p(x |x)= p(x |z,x)p(z|x)dz_E_z p(·|x)[p(x |z,x)]≈ p(x |z,x), (4) where the approximation is done by taking a single Monte Carlo sample from the expectation. Effectively, this means that only the probability over the relevant part of the output should be used, not the probability over the entire output. Intuitively, a mundane program generated following a very unlikely reasoning trajectory should have a reasonable generation probability, instead of being artificially down-weighted due to the unlikely reasoning that led to the program. 5.2 Length-Respecting Uniform Prior When deriving DME we initially assumed the batch’s samples are sampled from a uniform distribution over the space, hence the need for the 1/p(x′|x)1/p(x |x) importance sampling correction. However, which uniform distribution to use was not explicitly specified. Naïvely selecting a uniform distribution over sequences artificially up-weighs long sequences, e.g. for a vocabulary of size V there are V sequences with one token but V2V^2 with two. To reduce this length dependence for problems where shorter solutions are likely preferred, we define a uniform prior using the following process. Assume tokens are sampled uniformly from the vocabulary, and a sequence ends when an end-of-sequence token is sampled, and assume there is one such token in the vocabulary. Thus, the probability of sampling the token sequence making up x is U(x)=1V|x|U(x)= 1V^|x|, as each token in the sequence has probability 1/|V|1/|V| of being sampled. 5.3 Modified Importance Sampling Weights Altogether, these two modifications change the importance weights as follows. The path marginalization results in using w=R(x′)βp(x′|z,x)w= R(x )^βp(x |z,x) instead of w=R(x′)βp(x′|x)w= R(x )^βp(x |x), while the length-respecting uniform prior changes the numerator from w=R(x′)βp(x′|z,x)w= R(x )^βp(x |z,x) to w=R(x′)βU(x′)p(x′|z,x)=R(x′)βV|x′|p(x′|z,x)w= R(x )^βU(x )p(x |z,x)= R(x )^βV^|x |p(x |z,x). Note that the previous importance sampling weights could omit the U(x)U(x) term as it was constant. The path marginalization also yields a different heuristic interpretation of DME. When β→∞β→∞ and the generation probabilities between different samples are very different, DME simply implements a tie-breaker for hill climbing, where the highest reward sample with the lowest probability should be picked. The path marginalization results in the probability being calculated only over the relevant part of the answer instead of the entire generated sequence. Unless specified otherwise, the path marginalization is always used by default, while the default prior is the regular uniform distribution. 6 Experiments In this section we test DME’s sample efficiency on several problems. Specifically, we ask what is the probability DME and other methods can solve a problem given some number of samples. For difficult discovery problems requiring many samples any method which mode collapses or struggles sampling low probability solutions will have a low probability of success. 6.1 Regular Expressions To facilitate comparing DME to many existing methods we begin by testing its performance in a small-scale setting, generating regular expressions (regex). A regex is a short expression which accepts or rejects strings, akin to how a program can pass or fail unit tests. Regex are relatively high entropy strings, making them a good test case for rare event sampling.44 4 For example, a regex that accepts gmail addresses is ^[ .+\-]+@gmail\.com$ (Stelmaszczyk, 2013). A sampled regex can count as a success if it accepts a set of positive example strings (e.g. email addresses) and rejects a set of negative examples (e.g. malformed addresses). This setup is designed to be relatively simple and quick to run, enabling reliable probability estimates when comparing different methods. Specifically, we use Qwen3-0.6B as the generative model for all search methods and cap its output to a maximum of 16 tokens, due to the regex being relatively short. All methods are given 100k samples to find the regex and run on the same hardware – a single TitanRTX GPU – to facilitate runtime comparisons. The following methods are compared. Note that most of these methods perform local reward maximization and accordingly have objectives which permit mode collapse. First, we test three gradient-based discovery methods. Risk-sensitive GRPO (RS-GRPO) is a reinforcement-learning method with a diversity-inducing objective (Jiang et al., 2025). Test-time training (T-Discover) is a more elaborate RL-based discovery framework (Yuksekgonul et al., 2026). GFlowNets is a framework for optimizing similar distribution-matching objectives using methods from variational inference and reinforcement learning (Malkin et al., 2022b; Malkin et al., 2022a; Tiapkin et al., 2024). The remaining five baselines and DME are all gradient-free. In Best-of-N (BoN), also known as random sampling, the model is sampled IID while prompted to find the regex. Sequential conditioned sampling (SCS) is a strong baseline for discovery methods which uses a meta-model (also Qwen3-0.6B) to propose new prompts to find the regex (Gideoni et al., 2026). GEPA is a similar, more sophisticated method which proposes prompts also based on returned text-based task-specific information (Agrawal et al., 2026). (1+N) and (1,N) are two standard hill climbing evolutionary strategies. DME is run with β=1000β=1000. All methods had their setups and hyperparameters tuned, use a batch size of 64 where relevant, and where possible use the same prompts. Methods that rely on a reward signal use eaccuracye^accuracy as their reward function, where the accuracy is the percent of test cases which passed. Additional details are in Appendix B. For the test data we use the regex dataset of Siddiq et al. (2024) and filter to 20 regexes for which best-of-N gets between a 5% and 50% success rate, to ensure we focus on regexes which are difficult but solvable. All methods are run 10 times per regex except for some especially slow baselines, which are run only 3 times. We report the average success rate across regexes and 95% MOVER-aggregated confidence intervals (Zou and Donner, 2008). Figure 2: Average success rate across regexes for different methods. Black lines denote means, shaded regions are 95% confidence intervals. RS-GRPO T-Discover GFlowNets BoN SCS GEPA (1+N) (1,N) DME Time per run (minutes) 154.8±41.5154.8± 41.5 324.3±56.4324.3± 56.4 238.8±92.6238.8± 92.6 3.0±0.43.0± 0.4 18.1±2.918.1± 2.9 47.5±14.947.5± 14.9 6.7±1.96.7± 1.9 6.5±0.96.5± 0.9 7.2±0.57.2± 0.5 Markovian ✓ ✓ ✓ ✓ ✓ ✓ ✓ Requires a meta-model ✓ ✓ Requires weight updates ✓ ✓ ✓ Requires log-probs ✓ ✓ ✓ ✓ Table 1: Properties of different methods. DME is relatively fast, Markovian, does not need auxiliary models, and is inference only. However, it does require access to the generative model’s (potentially unnormalized) density. Runtime uncertainties are standard deviations. 6.1.1 Results Results across methods are shown in Figure 2. DME gets the highest success rate, almost twice the success rate of best-of-N. (1,N)(1,N) hill climbing is perhaps the only method other than DME to perform notably better than best-of-N, being effectively a simplified greedy heuristic version of DME in the β→∞β→∞ limit. We speculate (1,N)(1,N) hill climbing outperforms (1+N)(1+N) due to always refreshing the regex used as hint for generating the next batch, thus making its search more varied. (1+N)(1+N) changes its hint only if it generates an improvement over the best regex seen so far. It is surprising that GFlowNets have a low success rate given they are the only method other than DME which uses a distribution matching objective. We speculate GFlowNets’ low success is due to the gradient-based methods being unstable and mode collapsing even given distribution matching objectives. Empirically, the three tested baselines all work best only given very low learning rates – 10−710^-7 for T-Discover and GFlowNets, 10−610^-6 for RS-GRPO – with higher learning rates leading all three to mode collapse. 6.2 Coding Competitions For a more difficult, real-world task, we measure how well DME can solve problems from coding competitions. Specifically, we measure to what extent DME can enable a model to find correct answers which would otherwise require many samples to generate. We take a random 300 question subset of LiveCodeBench (Jain et al., 2025) and compare DME to (1,N)(1,N) and (1+N)(1+N) evolutionary strategies as they performed best for the regex task. DME uses β=1000β=1000 and the length-respecting prior, as coding competition problems often have succinct solutions – see Appendix C for ablations and details. All methods have a batch size of 8. The proposal model is Qwen2.5-Coder-1.5B Instruct (Hui et al., 2024). As in the regex case, the reward is eaccuracye^accuracy, where a given solution’s accuracy is the percent of unit tests that pass. Results are shown in Figure 3. DME is empirically more sample efficient than the evolutionary strategies, requiring fewer programs to solve the same fraction of problems. For example, DME reaches a ∼30% 30\% solve rate given 5050 programs, whereas both evolutionary strategies require generating ∼100 100 programs per question to reach the same performance. Figure 3: Different methods’ solve-rate for competitive coding problems as a function of the number of generated programs per question. Shaded regions denote the standard deviation over 5 seeds. 7 Math Bounds AlphaEvolve (Novikov et al., 2025) demonstrated how code generation can be used to discover new mathematical bounds. For certain problems, finding a bound is analogous to a black-box optimization problem, where a verifier gets as input a list of numbers and outputs a corresponding mathematical bound. A language model can then either generate code which directly constructs this list of numbers or finds it using various optimization methods. We test DME’s sample efficiency on 9 of the problems analyzed in Novikov et al. (2025), comparing it to two strong open-source baselines benchmarked by Gideoni et al. (2026). The first baseline is best-of-N (BoN), while the second is a more sequential version thereof called sequential-conditioned-sampling (SCS). Gideoni et al. (2026) found SCS to be a strong baseline when compared to open-source versions of AlphaEvolve, like ShinkaEvolve (Lange et al., 2026). All methods are given a sample budget of 5000 programs, sampling from GPT-OSS 20B (Agarwal et al., 2025). GPT-OSS’ knowledge cutoff is about a year before AlphaEvolve came out, so Novikov et al. (2025)’s solutions cannot be memorized. SCS hyperparameters are kept at their proposed defaults. DME uses β=106β=10^6 and to reduce runtime uses four chains or, as they are known in evolutionary algorithms parlance, four islands. All methods use the same prompt templates, and a temperature and top-p of 1. Following Gideoni et al. (2026), methods are considered tied if their results are within floating point errors as per numpy.isclose. Problem BoN SCS DME # problems DME ≥ method 9/9 9/9 Average rank (↓ ) 2.6 1.8 1.0 First autocorr. ineq. (↓ ) 1.524 1.517 1.510 Second autocorr. ineq. (↑ ) 0.8779 0.8784 0.8868 Uncertainty ineq. (↓ ) 0.3521 0.3521 0.3521 Erdős’ min. overlap (↓ ) 0.3815 0.3812 0.3809 Sums/diffs of sets (↑ ) 1.110 1.114 1.117 Max–min dist. ratio (↓ ) 12.8892 12.8892 12.8893 Heilbronn triangles (↑ ) 0.0339 0.0360 0.0365 Kissing number in 11D (↑ ) 348 566 582 Circle packing (↑ ) 2.62304 2.63109 2.63598 Table 2: Best discovered bounds by different methods when given 5000 samples. DME either finds the best bound or ties for first place for all problems. Arrows denote whether higher/lower is better. For all problems except circle packing results are truncated to the first different significant digit. Best results are bolded, second best are underlined. Table 2 shows the results. Relative to the two baselines, DME finds the best bound for all problems, either finding the best outright or tying for first place with a different method. To see how DME performs over a range of budgets we compute the probability of improvement between DME and the two baselines as a function of the number of generated programs. This comparison asks “when picking a random problem, what is the probability DME will outperform a given method?”, with ties counting as each method having an equal probability of winning (Agarwal et al., 2021). Figure 4 shows that for almost all sampling budgets DME is likely to improve over both baselines. Figure 4: DME’s probability of improvement over both baselines as a function of the number of generated programs. For most sample budgets DME is likely to outperform both baselines. For per problem breakdowns see Appendix D.3. 8 Related Work DME is an MCMC-based method to optimize distribution-matching objectives. These same objectives can be directly optimized using variational inference or reinforcement learning if weight updates are acceptable. To illustrate, given a target distribution π and a parameterized distribution pθp_θ, the minimum of KL(pθ||π)=pθ[logpθ−logπ]=pθ[−logπ]−H(pθ)KL(p_θ||π)=E_p_θ[ p_θ- π]=E_p_θ[- π]-H(p_θ) (5) is achieved when π=pθπ=p_θ, where H(p)H(p) denotes p’s entropy. Note that akin to DME, this objective does not require π to be normalized. When substituting the DME target distribution π(x)∝R(x)βπ(x) R(x)^β and taking the KL’s gradient shows how optimizing this objective is equivalent to a certain form of soft reinforcement learning. Explicitly, with the DME target distribution, KL(pθ||π)=−βpθ[logR(x)]+logZ−H(pθ)KL(p_θ||π)=- _p_θ[ R(x)]+ Z-H(p_θ). Rescaling by 1β 1β and setting R′(x)≔logR(x),α≔1βR (x) R(x),α 1β, one finds that: −∇θαKL(pθ||π)=pθ[R′(x)∇θlogpθ]+α∇θH(pθ),- _θ (p_θ||π)=E_p_θ[R (x) _θ p_θ]+α _θH(p_θ), (6) which is exactly a soft reinforcement learning policy gradient (Ziebart, 2010; Levine, 2018). There are various ways to optimize KL(pθ||π)KL(p_θ||π) or other distribution matching objectives. The most closely related method to DME is GFlowNets, which propose a variety of losses to train a generative model to sample proportionally to some reward (Bengio et al., 2021; Malkin et al., 2022a; Madan et al., 2023). Beyond requiring weight updates, GFlowNets differ from DME in that some formulations require only samples from the generative model and not densities, the sampling is amortized, and in some cases GFlowNets are known to be more sample efficient than MCMC based methods (Bengio et al., 2021). Some other rare event sampling methods, such as the cross-entropy method (Rubinstein, 1997), optimize KL(π||pθ)KL(π||p_θ). Surprisingly, the cross-entropy method was found to perform worse than other gradient-based methods for the regex experiments in section 6.1, leading to its exclusion. DME is also related to some symbolic code optimization methods. Superoptimization is the process of finding an optimal computer program for some task using exhaustive search methods (Massalin, 1987). Schkufza et al. (2013) propose stochastic superoptimization, a method which uses MCMC over a finite instruction set to more efficiently find an optimal program. DME is in a sense generative stochastic superoptimization, as it similarly uses MCMC to find optimal programs, while relying on generative models to increase the sample efficiency further. DME is related to various MCMC algorithms. Most closely related is iterated-sampling-importance-resampling (i-SIR), which similarly samples N candidates and then resamples a proposal from them weighted by some importance sampling weights (Tjelmeland, 2004; Andrieu et al., 2010); for an overview, see Laitinen and Vihola (2025). i-SIR was originally proposed to more efficiently calculate expectations relative to traditional Metropolis-Hastings MCMC, by using both rejected and accepted candidates (Tjelmeland, 2004). Similarly, upon completion DME can return a sample that was originally rejected if it has the highest seen reward. Many prior works using generative models have found discoveries using various heuristic evolutionary algorithms (Romera-Paredes et al., 2024; Novikov et al., 2025; OpenAI, 2026). These discoveries are likely enabled by a combination of using strong underlying generative models and generating many samples, thus sampling lower probability completions. Supporting this claim, Gideoni et al. (2026); Tsoukalas et al. (2026) find that simple search heuristics perform comparably well to more complicated evolutionary pipelines given the same generative model and sample budgets. Gideoni et al. (2026) further note that a problem’s formulation can make it more amenable to automated search, so a better problem formulation can make a method more sample efficient. 9 Discussion There are several straightforward extensions to DME. Incorporating other methods from either the MCMC or evolutionary algorithms literature could prove fruitful. If using gradients is permissible, it would be interesting using Hamiltonian Monte Carlo as a potentially much more efficient variant of DME. Another interesting approach would be controlling to what extent DME upweighs low probability samples by tempering the generative model’s probability in the importance weights. This could be easily achieved by setting p→pαp→ p^α for some 0<α≤10<α≤ 1. All this is left for future work. It would be interesting to test distribution matching methods further on cases with noisy, surrogate, or otherwise inexact rewards. Surrogate rewards are especially important when a real experiment is either expensive or takes a long time, in which case discovery methods can help proposing potential experiments. The expensive to evaluate, real reward might be maximized by a sample which seems suboptimal given a cheap proxy. In this case it is important to propose a diverse set of candidates which have a high but not only a maximal reward per the proxy. While local reward-maximization methods would reward-hack the proxy, distribution matching methods can be tuned to sample per various reward distributions. One limitation of DME is it requiring (unnormalized) probability densities. For some generative models finding the density is expensive, i.e. diffusion/flow matching models require at least an extra forward pass per step for calculating a generation’s density (Grathwohl et al., 2019). Additionally, not all closed-source language models return their log-probabilities. While likely possible adapting DME for these settings, doing so would require deriving a new transition kernel and be potentially nontrivial. More broadly, all training-free methods have some clear limitations. Any discovery which has a sufficiently low probability, likely due to being extremely dissimilar to the model’s training data, will still require a potentially unreasonable number of samples to find. Practically, training-free discovery methods are better thought of as ways to expand what a model can find given some number of samples, with new discoveries then being added to the model’s training data. A later finetuned or retrained generative model will then have an effectively larger set of reachable discoveries. Such an iterative procedure of expanding the set of reachable discoveries is sometimes described as how science progresses (Kauffman, 2000). References Agarwal et al. (2021) R. Agarwal, M. Schwarzer, P. S. Castro, A. C. Courville, and M. Bellemare Deep reinforcement learning at the edge of the statistical precipice. Advances in neural information processing systems 34, p. 29304–29320. Cited by: §7. Agarwal et al. (2025) S. Agarwal, L. Ahmad, J. Ai, S. Altman, A. Applebaum, E. Arbus, R. K. Arora, Y. Bai, B. Baker, H. Bao, et al. Gpt-oss-120b & gpt-oss-20b model card. arXiv preprint arXiv:2508.10925. Cited by: §7. Agrawal et al. (2026) L. A. Agrawal, S. Tan, D. Soylu, N. Ziems, R. Khare, K. Opsahl-Ong, A. Singhvi, H. Shandilya, M. J. Ryan, M. Jiang, et al. Gepa: reflective prompt evolution can outperform reinforcement learning. In International Conference on Learning Representations, Vol. 2026, p. 8479–8565. Cited by: §B.2, §6.1. Andrieu et al. (2010) C. Andrieu, A. Doucet, and R. Holenstein Particle markov chain monte carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72 (3), p. 269–342. Cited by: §3.2, §8. Atkinson (2008) K. E. Atkinson An introduction to numerical analysis. John wiley & sons. Cited by: Appendix A. Bengio et al. (2021) E. Bengio, M. Jain, M. Korablyov, D. Precup, and Y. Bengio Flow network based generative models for non-iterative diverse candidate generation. Advances in neural information processing systems 34, p. 27381–27394. Cited by: §1, §1, §2.2, §2, §8. Beyer and Schwefel (2002) H. Beyer and H. Schwefel Evolution strategies–a comprehensive introduction. Natural computing 1 (1), p. 3–52. Cited by: §3.2. Brooks et al. (2011) S. Brooks, A. Gelman, G. Jones, and X. Meng Handbook of markov chain monte carlo. CRC press. Cited by: §1. Eijkelboom et al. (2024) F. Eijkelboom, G. Bartosh, C. A. Naesseth, M. Welling, and J. van de Meent Variational flow matching for graph generation. Advances in Neural Information Processing Systems 37, p. 11735–11764. Cited by: §1. Gideoni et al. (2026) Y. Gideoni, S. Risi, and Y. Gal Simple baselines are competitive with code evolution. arXiv preprint arXiv:2602.16805. Cited by: §B.2, §D.1, §6.1, §7, §7, §8. Goodfellow et al. (2014) I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio Generative adversarial nets. Advances in neural information processing systems 27. Cited by: §1. Grathwohl et al. (2019) W. Grathwohl, R. T. Q. Chen, J. Bettencourt, and D. Duvenaud Scalable reversible generative models with free-form continuous dynamics. In International Conference on Learning Representations, External Links: Link Cited by: §9. Griggs et al. (2025) T. Griggs, S. Hegde, E. Tang, S. Liu, S. Cao, D. Li, C. Ruan, P. Moritz, K. Hakhamaneshi, R. Liaw, A. Malik, M. Zaharia, J. E. Gonzalez, and I. Stoica Evolving skyrl into a highly-modular rl framework. Note: Notion Blog Cited by: §B.2. Hastings (1970) W. K. Hastings Monte carlo sampling methods using markov chains and their applications. Cited by: §2.3. Ho et al. (2020) J. Ho, A. Jain, and P. Abbeel Denoising diffusion probabilistic models. Advances in neural information processing systems 33, p. 6840–6851. Cited by: §1. Hui et al. (2024) B. Hui, J. Yang, Z. Cui, J. Yang, D. Liu, L. Zhang, T. Liu, J. Zhang, B. Yu, K. Lu, et al. Qwen2.5-coder technical report. arXiv preprint arXiv:2409.12186. Cited by: §6.2. Jain et al. (2025) N. Jain, A. Gu, W. Li, F. Yan, T. Zhang, S. Wang, A. Solar-Lezama, K. Sen, and I. Stoica Livecodebench: holistic and contamination free evaluation of large language models for code. In International Conference on Learning Representations, Vol. 2025, p. 58791–58831. Cited by: §6.2. Jiang et al. (2025) Y. Jiang, J. Huang, Y. Yuan, X. Mao, Y. Yue, Q. Zhao, and L. Yan Risk-sensitive rl for alleviating exploration dilemmas in large language models. arXiv preprint arXiv:2509.24261. Cited by: §B.2, §6.1. Kauffman (2000) S. Kauffman Investigations. Oxford University Press New York. Cited by: §9. Kingma and Welling (2013) D. P. Kingma and M. Welling Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114. Cited by: §1. Kirkpatrick et al. (1983) S. Kirkpatrick, C. D. Gelatt Jr, and M. P. Vecchi Optimization by simulated annealing. Science 220 (4598), p. 671–680. Cited by: §2.1, §2.2. Lab (2026) T. M. Lab Tinker. External Links: Link Cited by: §B.2. Laitinen and Vihola (2025) P. Laitinen and M. Vihola Iterated sampling importance resampling with adaptive number of proposals. arXiv preprint arXiv:2512.00220. Cited by: §8. Lange et al. (2026) R. Lange, Y. Imajuku, and E. Cetin Shinkaevolve: towards open-ended and sample-efficient program evolution. In International Conference on Learning Representations, Vol. 2026, p. 74026–74078. Cited by: §7. Levine (2018) S. Levine Reinforcement learning and control as probabilistic inference: tutorial and review. arXiv preprint arXiv:1805.00909. Cited by: §8. Madan et al. (2023) K. Madan, J. Rector-Brooks, M. Korablyov, E. Bengio, M. Jain, A. C. Nica, T. Bosc, Y. Bengio, and N. Malkin Learning gflownets from partial episodes for improved convergence and stability. In International Conference on Machine Learning, p. 23467–23483. Cited by: §B.2, §8. Malkin et al. (2022a) N. Malkin, M. Jain, E. Bengio, C. Sun, and Y. Bengio Trajectory balance: improved credit assignment in gflownets. Advances in Neural Information Processing Systems 35, p. 5955–5967. Cited by: §B.2, §6.1, §8. Malkin et al. (2022b) N. Malkin, S. Lahlou, T. Deleu, X. Ji, E. Hu, K. Everett, D. Zhang, and Y. Bengio GFlowNets and variational inference. arXiv preprint arXiv:2210.00580. Cited by: §1, §2.2, §6.1. Massalin (1987) H. Massalin Superoptimizer: a look at the smallest program. ACM SIGARCH Computer Architecture News 15 (5), p. 122–126. Cited by: §8. Metropolis et al. (1953) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller Equation of state calculations by fast computing machines. The journal of chemical physics 21 (6), p. 1087–1092. Cited by: §2.3. Novikov et al. (2025) A. Novikov, N. Vũ, M. Eisenberger, E. Dupont, P. Huang, A. Z. Wagner, S. Shirobokov, B. Kozlovskii, F. J. Ruiz, A. Mehrabian, et al. Alphaevolve: a coding agent for scientific and algorithmic discovery. arXiv preprint arXiv:2506.13131. Cited by: §D.1, §1, §7, §7, §7, §8. OpenAI (2026) OpenAI Ten advances in mathematics and theoretical computer science. Technical report OpenAI. External Links: Link Cited by: §8. Romera-Paredes et al. (2024) B. Romera-Paredes, M. Barekatain, A. Novikov, M. Balog, M. P. Kumar, E. Dupont, F. J. Ruiz, J. S. Ellenberg, P. Wang, O. Fawzi, et al. Mathematical discoveries from program search with large language models. Nature 625 (7995), p. 468–475. Cited by: §8. Rosenthal (1970) H. P. Rosenthal On the subspaces of l p (p>2p>2) spanned by sequences of independent random variables. Israel Journal of Mathematics 8 (3), p. 273–303. Cited by: Appendix A. Rubin (1987) D. B. Rubin A noniterative sampling/importance resampling alternative to the data augmentation algorithm for creating a few imputations when fractions of missing information are modest: the sir algorithm. Journal of the American Statistical Association 82 (398), p. 543–546. Cited by: §3.2. Rubinstein (1997) R. Y. Rubinstein Optimization of computer simulation models with rare events. European Journal of Operational Research 99 (1), p. 89–112. Cited by: §8. Schkufza et al. (2013) E. Schkufza, R. Sharma, and A. Aiken Stochastic superoptimization. ACM SIGARCH Computer Architecture News 41 (1), p. 305–316. Cited by: §8. Siddiq et al. (2024) M. L. Siddiq, J. Zhang, L. Roney, and J. C. Santos Re (gex— dos) eval: evaluating generated regular expressions and their proneness to dos attacks. In Proceedings of the 2024 ACM/IEEE 44th international conference on software engineering: new ideas and emerging results, p. 52–56. Cited by: §B.1, Table 3, §6.1. Stelmaszczyk (2013) A. StelmaszczykAnswer to “regular expression – validate gmail addresses”(Website) Note: Stack Overflow answer, posted 2013-04-24 External Links: Link Cited by: footnote 4. Tiapkin et al. (2024) D. Tiapkin, N. Morozov, A. Naumov, and D. P. Vetrov Generative flow networks as entropy-regularized rl. In International Conference on Artificial Intelligence and Statistics, p. 4213–4221. Cited by: §6.1. Tjelmeland (2004) H. Tjelmeland Using all metropolis–hastings proposals to estimate mean values. Statistics 4. Cited by: §8. Tsoukalas et al. (2026) G. Tsoukalas, A. Kovsharov, S. Shirobokov, A. Surina, M. Firsching, G. Bérczi, F. J. Ruiz, A. Suggala, A. Z. Wagner, E. Wieser, et al. Advancing mathematics research with ai-driven formal proof search. arXiv preprint arXiv:2605.22763. Cited by: §8. Yang et al. (2025) A. Yang, A. Li, B. Yang, B. Zhang, B. Hui, B. Zheng, B. Yu, C. Gao, C. Huang, C. Lv, et al. Qwen3 technical report. arXiv preprint arXiv:2505.09388. Cited by: §C.2. Yuksekgonul et al. (2026) M. Yuksekgonul, D. Koceja, X. Li, F. Bianchi, J. McCaleb, X. Wang, J. Kautz, Y. Choi, J. Zou, C. Guestrin, and Y. Sun Learning to discover at test time. In Forty-third International Conference on Machine Learning, External Links: Link Cited by: §B.2, §6.1. Ziebart (2010) B. D. Ziebart Modeling purposeful adaptive behavior with the principle of maximum causal entropy. Carnegie Mellon University. Cited by: §8. Zou and Donner (2008) G. Y. Zou and A. Donner Construction of confidence limits about effect measures: a general approach. Statistics in medicine 27 (10), p. 1693–1702. Cited by: §6.1. Appendix A Proofs The formal version of theorem 1, including the necessary assumptions on the various distributions, is Theorem 1 (First order correction to q(x′|x)q(x |x)). Define the transition kernel q(x′|x)q(x |x) using the sampling-resampling process given in Algorithm 3. Assume that p(x′|x)p(x |x) has full support, that for all x the χ2χ^2-divergence χ2(π||p(⋅|x))=∫(π(x′)−p(x′|x)p(x′|x))2p(x′|x)dx′χ^2(π||p(·|x))= ( π(x )-p(x |x)p(x |x) )^2p(x |x)dx is finite, and that the importance weight’s third absolute moment is finite, so −i[|wj−Z|3]<∞E_-i[|w_j-Z|^3]<∞. Then, one has that q(x′|x)=π(x′)(1+D(x′|x)N)+O(N−3/2),q(x |x)=π(x ) (1+ D(x |x)N )+O (N^-3/2 ), (7) where D(x′|x)≔∫π(y)2p(y|x)y−π(x′)p(x′|x)D(x |x) π(y)^2p(y|x)dy- π(x )p(x |x). Proof of Theorem 1, First order correction to q(x′|x)q(x |x). We wish to find an exact analytic form for q(x′|x)q(x |x) based on the process used in DME, i.e. Algorithm 3. For simplicity we denote xix_i simply as x. Given x1′,…,xN′\x _1,...,x _N\, the probability of sampling x1′x _1 given each sample has weight wi≔R(xi′)βp(xi′|x)w_i R(x_i )^βp(x _i|x) is, as before, pi≔wi∑j=1Nwjp_i w_i _j=1^Nw_j. Thus, as xi′∼p(⋅|x)x _i p(·|x), the probability of sampling xi′x _i from the joint sampling process is Npip(xi′|x)Np_ip(x _i|x), where the N factor comes from the N samples being exchangeable. Note that pip_i implicitly still relies on the entire set x1,…,xN\x_1,...,x_N\. Throughout the proof we will take many expectations of the form x1,…,xi−1,xi+1,…,xN∼p(⋅|x)[⋅]E_x_1,...,x_i-1,x_i+1,...,x_N p(·|x)[·]. For convenience, denote −i≔x1,…,xi−1,xi+1,…,xN∼p(⋅|x)E_-i _x_1,...,x_i-1,x_i+1,...,x_N p(·|x). Denoting x′=xi′x =x _i and marginalizing over the remaining xj′x _j gives the exact x→x′x→ x transition probability: q(x′|x)=Np(x′|x)−i[pi]=NR(x′)β−i[1wi+∑j=1,j≠iNwj],q(x |x)=Np(x |x)E_-i[p_i]=NR(x )^βE_-i [ 1w_i+ _j=1,j =i^Nw_j ], (8) where the expectation is over all j≠ij =i. This form of the transition probability is difficult to use as the expectation is over a fraction. To reach an analytic formula we wish to express the fraction as an expansion, so the expectation can act on polynomial terms of wjw_j. As a prologue to such an expansion, note that: −i[wj]=∫R(xj′)βp(xj′|x)p(xj′|x)dxj′=∫R(xj′)βdxj′=Z,E_-i[w_j]= R(x _j)^βp(x _j|x)p(x _j|x)dx _j= R(x _j)^βdx _j=Z, (9) where recall that Z is the normalization constant of the target distribution, π(x)π(x). Thus, denote S≔(N−1)ZS (N-1)Z and Δ≔∑j=1,j≠iNwj−S _j=1,j =i^Nw_j-S. Note that x′j∼p(⋅|x),j≠i[Δ]=0E_x _j p(·|x),j =i[ ]=0, by design. The rest of the proof boils down to straightforward if somewhat tedious algebra and a truncated Taylor expansion around Δwi+S=0 w_i+S=0. Using the new notation, one has that: −i[1wi+∑j=1,j≠iNwj]=−i[1wi+S+Δ]=1wi+S−i[11+Δwi+S].E_-i [ 1w_i+ _j=1,j =i^Nw_j ]=E_-i [ 1w_i+S+ ]= 1w_i+SE_-i [ 11+ w_i+S ]. (10) Now, note that as wi+S+Δ=∑j=1Nwj>0w_i+S+ = _j=1^Nw_j>0 then the denominator is always positive and the expectation has no poles. Taking an exact finite Taylor expansion around Δwi+S=0 w_i+S=0 gives: (1+Δwi+S)−1=1−Δwi+S+(Δwi+S)2−(Δwi+S)31+Δwi+S, (1+ w_i+S )^-1=1- w_i+S+ ( w_i+S )^2- ( w_i+S )^31+ w_i+S, (11) where we used the expansion (1+x)−1=1−x+x2−x31+x(1+x)^-1=1-x+x^2- x^31+x [Atkinson, 2008, Eq. (1.1.8)]. Plugging the first two terms back into the expectation in equation 10 gives: −i[1]=1,−i[Δwi+S]=0,E_-i [1 ]=1, _-i [ w_i+S ]=0, (12) as −i[Δ]=0E_-i[ ]=0. For the third term, we have that: −i[Δ2]=−i[∑j,k=1,j,k≠iNwjwk]−S2.E_-i [ ^2 ]=E_-i [ _j,k=1,j,k =i^Nw_jw_k ]-S^2. (13) For the (N−1)2−(N−1)(N-1)^2-(N-1) terms where j≠kj =k one has that −i[wjwk]=−i[wj]−i[wk]=Z2,E_-i [w_jw_k ]=E_-i [w_j ]E_-i [w_k ]=Z^2, which leaves only the N−1N-1 remaining −i[wj2]E_-i [w_j^2 ] terms. Note that: −i[wj2]=∫R(xj′)2βp(xj′|x)2p(xj′|x)dxj′=Z2∫π(y)2p(y|x)y,E_-i [w_j^2 ]= R(x _j)^2βp(x _j|x)^2p(x _j|x)dx _j=Z^2 π(y)^2p(y|x)dy, (14) where for convenience the integrand is denoted as y instead of xj′x _j. Plugging these results back into equation 13 gives: −i[Δ2] _-i [ ^2 ] =((N−1)2−(N−1))Z2+(N−1)−i[wj2]−S2 = ((N-1)^2-(N-1) )Z^2+(N-1)E_-i [w_j^2 ]-S^2 (15) =(N−1)Z2(N−2+∫π(y)2p(y|x)dy−(N−1)) =(N-1)Z^2 (N-2+ π(y)^2p(y|x)\,dy-(N-1) ) =(N−1)Z2(∫π(y)2p(y|x)dy−1) =(N-1)Z^2 ( π(y)^2p(y|x)\,dy-1 ) Now, temporarily denoting I≔∫π(y)2p(y|x)y−1I π(y)^2p(y|x)\,dy-1, one has that the expectation of the second order term in equation 11 is: −i[(Δwi+S)2]=(N−1)Z2I(wi+S)2.E_-i [ ( w_i+S )^2 ]= (N-1)Z^2I (w_i+S )^2. (16) Specifically, note that χ2(π||p(⋅|x))=Iχ^2(π||p(·|x))=I, which is assumed to be finite and hence well defined. In the large N limit it is clear that (N−1)Z2(wi+S)2≈1N (N-1)Z^2(w_i+S)^2≈ 1N. Thus, this expression can be slightly cleaned up if discarding higher order terms in N. More precisely, note that: (N−1)Z2(wi+S)2=(N−1)Z2(wi+(N−1)Z)2=1N+O(N−2), (N-1)Z^2(w_i+S)^2= (N-1)Z^2(w_i+(N-1)Z)^2= 1N+O(N^-2), (17) as wi=O(1)w_i=O(1). Thus, equation 16 can be simplified to: −i[(Δwi+S)2]=IN+O(N−2).E_-i [ ( w_i+S )^2 ]= IN+O(N^-2). (18) Now all that remains is to bound the third order remainder term in equation 11 and aggregate results. Specifically, we wish to bound: |−i[(Δwi+S)31+Δwi+S]|<−i[|(Δwi+S)31+Δwi+S|]=−i[|Δ3(wi+S)2(wi+S+Δ)|]<−i[|Δ3S2(wi+S+Δ)|]. |E_-i [ ( w_i+S )^31+ w_i+S ] |<E_-i [ | ( w_i+S )^31+ w_i+S | ]=E_-i [ | ^3(w_i+S)^2(w_i+S+ ) | ]<E_-i [ | ^3S^2(w_i+S+ ) | ]. (19) This last expression is nontrivial to bound as wi+S+Δw_i+S+ , while necessarily positive, could be close to zero. To bound it we split the expectation into cases, where either Δ≥−S/2 ≥-S/2 or Δ<−S/2 <-S/2. The first case is “good” in that the denominator is clearly nonzero, so we shall denote it with the indicator function G 1_G. Similarly, the second case is“bad” in that the denominator could be close to 00, so we will denote it with B 1_B, where clearly G+B=1 1_G+ 1_B=1 as the two cases are complementary. In the good case we have that: |Δ3S2(wi+S+Δ)|≤2|ΔS|3. | ^3S^2(w_i+S+ ) |≤ 2 | S |^3. (20) Now, we wish to bound |ΔS|3 | S |^3. As the wjw_j are IID and we assume finite third absolute moments (−i[|wj−Z|3]<∞E_-i[|w_j-Z|^3]<∞), for some constant C Rosenthal’s inequality [Rosenthal, 1970] gives: −i[|Δ|3]≤C(∑j≠i−i[|wj−Z|3]+(∑j≠i−i[(wj−Z)2])3/2),E_-i[| |^3]≤ C ( _j =iE_-i[|w_j-Z|^3]+ ( _j =iE_-i[(w_j-Z)^2] )^3/2 ), (21) where note that each sum has O(N)O(N) terms. Thus, ignoring constant factors: −i[|Δ|3]≤C((N−1)O(1)+((N−1)O(1))3/2)=C(O(N)+(O(N))3/2)=O(N3/2),E_-i[| |^3]≤ C ((N-1)O(1)+ ((N-1)O(1) )^3/2 )=C (O(N)+ (O(N) )^3/2 )=O(N^3/2), (22) which as S=(N−1)ZS=(N-1)Z implies −i[|Δ/S|3]=O(N−3/2)E_-i[| /S|^3]=O(N^-3/2) for the part of the expectation where Δ≥−S/2 ≥-S/2. As for the bad region, we shall show that the probability of Δ<−S/2 <-S/2 is very small, so the contribution from the bad region is negligible. Explicitly, note that −S/2>Δ=∑j≠iwj−S>−S-S/2> = _j =iw_j-S>-S as wj>0w_j>0. Thus, |Δ|<S| |<S and we have that: |Δ3S2(wi+S+Δ)|<S3S2wi=Swi=O(N), | ^3S^2(w_i+S+ ) |< S^3S^2w_i= Sw_i=O(N), (23) as wiw_i is constant and S=(N−1)Z=O(N)S=(N-1)Z=O(N). Thus, we have that −i[B|Δ3S2(wi+S+Δ)|]<O(N)−i[B]=O(N)Pr(Δ<−S/2)E_-i [ 1_B | ^3S^2(w_i+S+ ) | ]<O(N)E_-i[ 1_B]=O(N)Pr( <-S/2) (24) Thus, we wish to find a sufficiently tight bound for Pr(Δ<−S/2)Pr( <-S/2). Note that Δ<−S/2⟹|Δ|>S/2 <-S/2 | |>S/2. Intuitively, for large N the probability of the sum of displacements (Δ ) being larger than its mean by a constant factor (in this case 12 12) is extremely low. Practically, we bound Pr(Δ<−S/2)Pr( <-S/2) using Hoeffding’s inequality. Hoeffding’s inequality requires bounded random variables, but wjw_j is generically unbounded. Thus, for some M>0M>0 define yj≔min(wj,M)y_j (w_j,M). As wj>0w_j>0 we have that 0<yj≤M0<y_j≤ M and is thus bounded. As by definition yj≤wjy_j≤ w_j we have that Δ+S=∑j≠iwj<S/2⇒∑j≠iyj<S/2 +S= _j =iw_j<S/2 _j =iy_j<S/2. Moreover, note that Pr(∑j≠iwj<S/2)≤Pr(∑j≠iyj<S/2)Pr( _j =iw_j<S/2) ( _j =iy_j<S/2), as if the first inequality holds (∑j≠iwj<S/2 _j =iw_j<S/2) the second inequality (∑j≠iyj<S/2 _j =iy_j<S/2) holds as well. Recall now that M is arbitrary, thus we can define M to be sufficiently large such that −i[yj]≥34ZE_-i[y_j]≥ 34Z. The 34 34 is arbitrary; any fraction larger than 12 12 would suffice, as Hoeffding’s inequality bounds the probability that a random variable would deviate from its mean by a constant factor. Thus, applying Hoeffding’s inequality on the sum of yjy_j implies that: Pr(Δ<−S/2)=Pr(∑j≠iwj<S/2)≤Pr(∑j≠iyj<S/2)≤e−cNPr( <-S/2)=Pr ( _j =iw_j<S/2 ) ( _j =iy_j<S/2 )≤ e^-cN (25) for some c>0c>0. Thus, we have that: −i[B|Δ3S2(wi+S+Δ)|]<O(N)Pr(Δ<−S/2)≤O(N)e−cN=o(N−3/2).E_-i [ 1_B | ^3S^2(w_i+S+ ) | ]<O(N)Pr( <-S/2)≤ O(N)e^-cN=o(N^-3/2). (26) Technically this expression equals o(N−k)o(N^-k) for any positive k, but 3/23/2 suffices here. Altogether, the results for the good and the bad regions imply that: |−i[(Δwi+S)31+Δwi+S]| |E_-i [ ( w_i+S )^31+ w_i+S ] | <−i[|Δ3S2(wi+S+Δ)|] <E_-i [ | ^3S^2(w_i+S+ ) | ] (27) <2−i[G|Δ/S|3]+−i[B|Δ3S2(wi+S+Δ)|] <2E_-i [ 1_G| /S|^3 ]+E_-i [ 1_B | ^3S^2(w_i+S+ ) | ] <2−i[|Δ/S|3]+o(N−3/2) <2E_-i [| /S|^3 ]+o(N^-3/2) =O(N−3/2)+o(N−3/2) =O(N^-3/2)+o(N^-3/2) =O(N−3/2). =O(N^-3/2). thus bounding the remainder term. Taking the results of equations 18 and 27, plugging them into equation 11 and taking an expectation gives: −i[(1+Δwi+S)−1]=1+IN+O(N−3/2).E_-i [ (1+ w_i+S )^-1 ]=1+ IN+O(N^-3/2). (28) Plugging results back in now until we reach the original expression for q(x′|x)q(x |x), plugging equation 28 into equation 10 gives: −i[1wi+∑j=1,j≠iNwj]=1wi+S(1+IN+O(N−3/2)).E_-i [ 1w_i+ _j=1,j =i^Nw_j ]= 1w_i+S (1+ IN+O(N^-3/2) ). (29) Now, note that: 1wi+S=1NZ11+wi−ZNZ=1NZ11−1−wi/ZN=1NZ(1+1−wi/ZN+O(N−2)) 1w_i+S= 1NZ 11+ w_i-ZNZ= 1NZ 11- 1-w_i/ZN= 1NZ (1+ 1-w_i/ZN+O(N^-2) ) (30) Thus, one has that: −i[1wi+∑j=1,j≠iNwj] _-i [ 1w_i+ _j=1,j =i^Nw_j ] =1NZ(1+O(N−3/2)+IN+1−wi/ZN+O(N−2)) = 1NZ (1+O(N^-3/2)+ IN+ 1-w_i/ZN+O(N^-2) ) (31) =1NZ(1+I+1−π(x′)p(x′|x)N+O(N−3/2)), = 1NZ (1+ I+1- π(x )p(x |x)N+O(N^-3/2) ), where the O(N−2)O(N^-2) term is absorbed by the O(N−3/2)O(N^-3/2). Plugging 31 back into equation 8 gives: q(x′|x) q(x |x) =NR(x′)β−i[1wi+∑j=1j≠iNwj] =NR(x )^βE_-i [ 1w_i+ _ subarraycj=1\\ j≠ i subarray^Nw_j ] (32) =π(x′)(1+I+1−π(x′)p(x′|x)N)+O(N−3/2) =π(x ) (1+ I+1- π(x )p(x |x)N )+O(N^-3/2) =π(x′)(1+1N(∫π(y)2p(y|x)dy−π(x′)p(x′|x)))+O(N−3/2). =π(x ) (1+ 1N ( π(y)^2p(y|x)dy- π(x )p(x |x) ) )+O(N^-3/2). Defining D(x′|x)≔∫π(y)2p(y|x)y−π(x′)p(x′|x)D(x |x) π(y)^2p(y|x)dy- π(x )p(x |x) thus finishes the proof. ∎ The proof of corollary 2 is relatively straightforward. As the version stated in the main body is precise, we give only the proof here. Proof of Corollary 2, First order correction as a χ2χ^2 divergence. Directly calculating: x′∼p(⋅|x)[D(x′|x)] _x p(·|x)[D(x |x)] =∫π(y)2p(y|x)−π(y)p(y|x)p(y|x)y = π(y)^2p(y|x)- π(y)p(y|x)p(y|x)\,dy (33) =∫π(y)2p(y|x)dy−1 = π(y)^2p(y|x)\,dy-1 =∫π(y)2p(y|x)−p(y|x)y = π(y)^2p(y|x)-p(y|x)\,dy =∫(π(y)−p(y|x)p(y|x))2p(y|x)dy. = ( π(y)-p(y|x)p(y|x) )^2p(y|x)\,dy. where we used that ∫π(y)y=∫p(y|x)y=1 π(y)dy= p(y|x)dy=1. This equalling χ2(π||p(⋅|x))χ^2(π||p(·|x)) is by definition of the χ2χ^2 divergence. ∎ Appendix B Regex Details B.1 Example Prompts Each regex in the dataset of Siddiq et al. [2024] has a prompt which briefly explains what the regex does and gives some examples for matches and non-matches. Table 3 contains these prompts and the dataset’s reference regexes for the three regexes used as a validation set. Note that the prompts only show a few of the examples, while the accuracy is measured over all 25 (13 matches and 12 non-matches). All methods which get as input some prior regex as a hint used the prompt Please output a regular expression that detects the following, with nothing preceding or succeeding it. description The last proposed regex was previous regex. The regex is: where description is the regex’ prompt from the dataset and previous regex is the hint, which is empty at the first step. Methods with no hint omit the sentence with The last proposed regex was. Example correct regex Description +([^.,])?$ Accepts positive integers. Match examples: - “34” - “1” - “1000000000” Non-match examples: - “-34” - “-1” +\/? *$ Used to validate fractions or UK odds. It will match integers, allowing 10/1 to be entered as 10, and fractions. It will not match negative numbers or fractions, spaces or any alphabetical characters. Match examples: - “10” - “2/3” - “100/30” Non-match examples: - “-10” - “-2/3” ^-? +(\. +)?$ Matches any unsigned or signed floating point number/numeric string. Match examples: - “123” - “3.14159” - “-1.2” Non-match examples: - “abc” - “3.4.2” Table 3: Prompts (refined_prompt in the dataset) and reference regexes for the three validation regexes. References are the ground-truths given in the dataset of Siddiq et al. [2024]. B.2 Regex Baselines In this section we review the different baselines used for the regex discovery task and how their setups and/or hyperparameters were tuned. Baselines not mentioned here do not have any hyperparameters and use the same prompts as DME. For all methods we used 3 regexes not included in the test set as validation regexes and used the hyperparameters that had the highest average success rate over these 3. For RS-GRPO [Jiang et al., 2025] we sweeped their β parameter and the learning rate. Note that the RS-GRPO objective is of the form 1βlog[eβR(x)] 1β [e^β R(x)], so their β is different to the one used by other methods here. We ran a hyperparameter sweep over β=2,8β=2,8, as these are the two β used for the main results in Jiang et al. [2025], and learning rate=10−4,10−5,10−6=10^-4,10^-5,10^-6, finding that β=8,lr=10−6β=8,lr=10^-6 is best. For T-Discover [Yuksekgonul et al., 2026] we initially tried to use an open-source replication that does not require the Tinker API [Lab, 2026] but did not manage to get it to perform well on the regexes. Thus, we used the official T-Discover codebase in conjunction with skyRL tx [Griggs et al., 2025] to simulate the Tinker backend. We ran a hyperparameter sweep over learning rates 10−5,10−6,10−710^-5,10^-6,10^-7 and KL coefficients 0,0.1,0.010,0.1,0.01, finding that 10−710^-7 and 00 respectively performed best. For GFlowNets we implemented versions using both the trajectory balance loss [Malkin et al., 2022a] and the subtrajectory balance loss [Madan et al., 2023], finding that the latter gave slightly better performance. The subtrajectory balance loss slightly increases runtime relative to the trajectory balance loss as the former requires evaluating the reward given to all partial sequences and not only at the terminal state, which in our case is the final completion. We ran a hyperparameter sweep over β=0,1,3,10,30,100,300,1000,3000β=0,1,3,10,30,100,300,1000,3000 and learning rates 10−6,10−7,10−810^-6,10^-7,10^-8, finding that 1000,10−71000,10^-7 respectively performed best, although many β and learning rates performed equally well. For GEPA [Agrawal et al., 2026] we manually iterated over the prompt and additional textual information given to the model. Moreover, we force the model to stop thinking and give an answer after 40964096 tokens as otherwise the model tends to exhaust its context length without outputting a new prompt. Otherwise, GEPA had no hyperparameters tuned. SCS [Gideoni et al., 2026] used the same prompt as GEPA, only a single trial, and otherwise uses the same hyperparameters as given by Gideoni et al. [2026]. B.3 Ablations In this section we test how different design choices, namely DME’s beta parameter and the two tweaks to the base algorithm, affect its performance on the regex task. All experiments and reported performance are over the same 20 test regexes. To test how DME’s performance is affected by beta we sweep over betas from 10010^0 to 10610^6. Figure 5 shows that above about β=103β=10^3 DME on average performs the same. Note that this does not mean that DME’s performance monotonically increases for all problems, e.g. for the third validation regex in Table 3 β=100β=100 has a success rate of 96% while the success rate of β=3000β=3000 is 86%. Figure 5: DME’s average performance for the regex task increases with β, stabilizing after β=103β=10^3. Uncertainties are 95% confidence intervals. For the regex task it does not seem to matter much whether the path marginalization or length-respecting prior are used, see Table 4. This makes sense as the output is relatively short and does not contain long, explicit reasoning leading up to it. Path Marginalisation ✓ ✗ Length-Respecting Prior ✓ 45.5% 49.5% ✗ 46.5% 48.5% Table 4: DME average success rate (β=1000β=1000), with and without the length-respecting prior and path marginalisation. 95% confidence intervals are about ±5%± 5\% for all cells. Appendix C Coding Competitions Details C.1 DME Ablations and Hyperparameter Selection To initially select which hyperparameters to use for DME we ran it on a random subset of 100 LiveCodeBench problems, distinct and harder than those used for comparing methods in Figure 3. We sweep over β=101,103,106β=10^1,10^3,10^6 and test DME with vs without the length-respecting uniform prior. Figure 6 shows the result of that sweep, where initially promising runs were run with an additional seed. β=103β=10^3 seems to work better than β=101β=10^1 or β=106β=10^6, albeit only when using the length-respecting uniform prior. Figure 6: DME for different betas and with different priors, over a hard 100 question subset of LiveCodeBench. Configurations with shaded regions report mean and standard deviation over two seeds. C.2 Results for a Reasoning Model We compare DME to (1,N)(1,N) and (1+N)(1+N) when using a reasoning model, namely Qwen3-1.7B [Yang et al., 2025]. DME here uses β=106β=10^6. Surprisingly, in this case we find DME to slightly underperform both evolutionary strategies. As the results for section 7 show, DME can be competitive when applied to reasoning models over many different problems, so it is unclear why DME underperformed here. Figure 7: Different methods’ performance when using a reasoning model as the proposer. Surprisingly, here DME slightly underperforms both evolutionary strategies. Appendix D Math Bounds Details D.1 DME Details and Hyperparameter Tuning DME’s reward is ebounde^bound, with the bound’s sign flipped when it needs to be minimized rather than maximized. Note that both BoN and SCS are invariant to the reward parameterization as they have no explicit fitness-based selection during the search, only selecting the min/max discovered bound after evaluating all programs. DME’s hyperparameters were tuned by running short, single-chain runs that generate up to 1000 programs on the third autocorrelation inequality problem from the original AlphaEvolve paper [Novikov et al., 2025]. This problem was not studied in Gideoni et al. [2026], hence it being excluded from the main results. Table 5 shows the results found per run. Both using the standard prior and increasing beta seemed to help. Length-respecting prior β ✓ ✗ 10310^3 1.556 1.508 10610^6 1.489 1.488 Table 5: DME hyperparameter tuning on the third autocorrelation inequality problem, lower is better. Each entry is the best discovered upper bound. D.2 Probability of Improvement Calculations Additional programs were generated to get a better estimate for the probability of improvement, with the average then being calculated by bootstrapping. Specifically, we sample an additional 1000 programs from both BoN and SCS. More samples are generated from SCS by generating more trials. We sample 3 additional chains for DME, each running for the same number of steps as the original 4 did. The bootstrapping consists of sampling from each method without replacement, namely from the BoN programs, the SCS trials, and the DME chains. D.3 Per-Problem Breakdown See Figure 8. Note that per-problem estimates are much less reliable than the aggregate estimate due to having less data. Figure 8: Per-problem DME probability of improvement.