LLM-Guided Evolutionary Program Synthesis for Quasi-Monte Carlo Design
Authors: Amir Sadikov
Venue: ICLR 2026 Poster (2026)
Keywords: Evolutionary Computation
Large Language Models (LLMs)
Quasi-Monte Carlo (QMC)
Star Discrepancy
Sobol' Sequences
Algorithmic Discovery
Quantitative Finance
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: LLM-guided evolutionary program synthesis discovers new best-known 2D finite low-discrepancy sets (N≥40), improves 3D benchmarks beyond proven optimality, and evolves Sobol’ parameters that reduce rQMC error in 32-D option pricing.

Abstract: Low-discrepancy point sets and digital sequences underpin quasi-Monte Carlo (QMC) methods for high-dimensional integration. We cast two long-standing QMC design problems as program synthesis and solve them with an LLM-guided evolutionary loop that mutates and selects code under task-specific fitness: (i) constructing finite 2D/3D point sets with low star discrepancy, and (ii) choosing Sobol’ direction numbers that minimize randomized quasi-Monte Carlo (rQMC) error on downstream integrands. Our two-phase procedure combines constructive code proposals with iterative numerical refinement. On finite sets, we rediscover known optima in small 2D cases and set new best-known 2D benchmarks for N ≥ 40, while matching known 3D optima up to the proven frontier (N ≤ 8) and reporting improved 3D benchmarks beyond. On digital sequences, evolving Sobol' parameters yields consistent reductions in rQMC mean-squared error for several 32-dimensional option-pricing tasks relative to widely used Joe–Kuo parameters, while preserving extensibility to any sample size and compatibility with standard randomizations. Taken together, the results demonstrate that LLM-driven evolutionary program synthesis can automate the discovery of high-quality QMC constructions, recovering classical designs where they are optimal and improving them where finite-N structure matters. Data and code are available at https://github.com/hockeyguy123/openevolve-star-discrepancy.

Published as a conference paper at ICLR 2026
LLM-GUIDED EVOLUTIONARY PROGRAM SYNTHESIS FOR QUASI-MONTE CARLO DESIGN
Amir Sadikov
University of California, San Francisco amir. sadikoveucsf.ed.
ABSTRACT
Low-discrepancy point sets and digital sequences underpin quasi-Monte Carlo (QMC) methods for high-dimensional integration. We cast two long-standing QMC design problems program synthesis and solve them with an LLM-guided evolutionary loop that mutates and selects code under task-specific fitness: (i) constructing finite 2D/3D point sets with low star discrepancy, and (ii) choos- ing Sobol' direction numbers that minimize randomized QMC error on down- stream integrands. Our two-phase procedure combines constructive code propos- als with iterative numerical refinement. On finite sets, we rediscover known op- tima in small 2D cases and set new best-known 2D benchmarks for N Z 40, while matching most known 3D optima up to the proven frontier (N S 8) and reporting improved 3D benchmarks beyond. On digital sequences, evolving Sobol' parameters yields consistent reductions in randomized quasi-Monte Carlo (rQMC) mean-squared error for several 32-dimensional option-pricing tasks rel- ative to widely used Joe-Kuo parameters, while preserving extensibility to any sample size and compatibility with standard randomizations. Taken together, the results demonstrate that LLM-driven evolutionary program synthesis auto- mate the discovery of high-quality QMC constructions, recovering classical de- signs where they are optimal and improving them where finite-N structure matters. Data and code are available at hhtpps: //tgthtttbtttt///// com coml//hockkyyyy123// openevolve star di scrrepancy
INTRODUCTION
Numerical integration in high dimensions is cornerstone of modern science and engineering. While standard Monte Carlo (MC) methods offer robust approach, their convergence rate, governed by the Central Limit Theorem, is often insufficient for applications requiring high precision (Glasserman,
2003). Quasi-Monte Carlo (QMC) methods provide compelling alternative by replacing pseu- dorandom samples with deterministic, highly uniform point sets (Dick & Pillichshammer [2010] Owen 1995). The uniformity of these is quantified by their discrepancy, with lower values corresponding to more evenly distributed points. The -Hlawka inequality provides the the- oretical underpinning for QMC, guaranteeing that the integration error is bounded by the product of the variation of the integrand and the star discrepancy of the point set (Koksma) 1964 Hlawka
I961].
This relationship has fueled decades of research into the discovery and construction of low- discrepancy sets. primary challenge in this field is the "ab initio" construction of finite point set of N points in d dimensions that minimizes star discrepancy. This is a problem of combinatorial complexity, and while optimal solutions have been found in 2D for N < 21 and in 3D for N s (Clément al.) 2024), the problem remains largely open for larger N and higher dimensions.
related challenge is the construction of infinite low-discrepancy sequences, such as those by Hal- ton, Hammersley, and Sobol'. Among these, Sobol' sequences are particularly prominent due to their excellent distribution properties and efficient generation (Sobol" [967]. However, their quality highly dependent on set of integer parameters known as direction numbers, and finding optimal parameters that ensure uniformity across all owwddimeniional projections a difficult combinatorial search problem (Joe& & Kuo 2008).

---

Published as a conference paper at ICLR 2026
Recent breakthroughs in Large Language Models (LLMs) have demonstrated their remarkable pabilities in code generation, logical reasoning, and pattern recognition (Gemini Team et al. (025). This has motivated the development of automated scientific discovery systems, such as AlphaE- volve, which leverage L.MSs to navigate complex search spaces (Novikov et al. 2025). In this work, utilize the OpenEvolve framework (Sharma 2025), open-source implementation based the principles of AlphaEvolve, to tackle the aforementioned challenges in discrepancy theory. We treat the construction of these mathematical objects as a program synthesis problem within an evo- lutionary framework (Koza 1994). The LLM acts an intelligent mutation operator, iteratively modifying code that generates candidate solutions based on feedback from a fitmess function.
This paper presents the successful application of this LLM-driven evolutionary approach to two fundamental problems:
1. Discovering finite 2D and 3D point sets with state-of-the-ee low star discrepancy.
2. Searching for superior Sobol' direction numbers to minimize randomized quasi-Mont Carlo (rQMC) integration error in the high-dimensional context of pricing exotic financial options (Paskov & Traub 1995).
Our results show that this methodology match and, in several important cases, surpass the best known human-derived solutions, without any task-specific training of the LLM. This suggests that LLM--riven evolutionary search is promising new paradigm for exploration and discovery in com- putational mathematics.
THEORETICAL BACKGROUND
2.1 STAR DISCREPANCY
Star discrepancy is the most common measure for quantifying the uniformity of point within the d-dimensional unit hypercube, [0, 1]dd (Owen, 1995). captures the largest deviation between the volume of axis-aligned "anchor box" and the fraction of points contained within it.
Definition 1 (Star Discrepancy). Let P =
x be a set of N points in [0, 1j2. An anchor box (0.q) for any q (91.... ....... qd) E [0, 1] is the hyperrectangle [0, q1) x [0, qa). The star discrepancy Dřy of the set P is defined as:
(PP))..)) Vol((o,
(P)
sup q
Here, the supremum is taken over all possible anchor boxes. A small D value implies that for any anchor box, the fraction of points falling within it is good approximation of its volume, indicating high uniformity. The practical importance of discrepancy for QMC is cemented by the Koksma- Hlawka inequality (Koksma 1964 Hlawka 1961). which bounds the error of QMC integration:
U f(u) du
f(x)) S VV() (P)
where V(J) is the total variation of the function f in the sense of Hardy and Krause (Niederreiter) 1992 Dick Pillichhaammer, 2010). This inequality guarantees that point sets with lower star discrepancy lead to smaller integration error bounds.
2.2 SOBOL' SEQUENCES AND DIRECTION NUMBERS
Sobol' sequences are class of low-discrepancy sequences that are particularly effective for QMC integration (Sobol") 1967 Dick & Pillichshammer 2010). They are constructed using the properties of primitive polynomials over the finite field two elements, Fg.
Definition 2 (Sobol' Sequence Construction). For each dimension j > 1, primitive polynomial over F2 of degree s, is chosen (Dick & Pillichshammer 2010
(2) :: + 01,j2"/-
152 +1
(3)
2

---

Published as a conference paper at ICLR 2026
where the coefficients dk.j either 0 1. From this polynomial, sequence of positive, odd integers called direction numbers ma.j (for lo = (5). are chosen freely.
For k Sj. subsequent direction numbers (expressed as fixed-point binary fractions Uk.j
mk. (2*) are generated by the canonical Sobol' recurrence (Joe & Kuo, 2008 Dick & Pillichsham- mer 2010
Uk.j 011 + 2k-25
1ky-115 (Uk-eyo > 8,)
(4)
where denotes bitwise XOR the binary fixed-poine representation and >> S is right bit-shift by 8y places. For j = 1, we set 2k 2 (van der Corput base-2).
The j-th coordinate of the i-th point in the sequence, ziji is then generated using Gray-code bits:
Ji.j Il Ukjj
with
122 g = i (i
Σ gk 2k2k
(5)
A21
The quality of the Sobol' sequence, particularly the uniformity of its low-dimensional projections, critically dependent on the choice of the primitive polynomials and the initial direction numbers (m1,.... m,). The work of Joe & Kuo (2008) provides widely used set of these parameters that serve strong baseline. In tables we encode the polynomial coefficients as,j as an integer Aj (binary bitmask).
3 RELATED WORK
Classical approaches for generating low-diserypancy sets are primarily number-theoretic. Founda- tional methods include Halton, Hammersley, and Sobol' sequences, which are designed achieve superior asymptotic uniformity compared to random sampling. While powerful, these classical con- structions are not always optimal for finite number of points N. Mathematical programming has been used to find provably optimal sets, though these approaches are computationally intensive and limited small instances [Clément et 2024). Heuristic methods, such as genetic algorithms and threshold accepting, have been applied to tackle larger instances by searching the space of point configurations (Clément al. 2023)..
More recently, machine learning techniques have been introduced to this domain. Message-Passing Monte Carlo (MPMC) leverages Graph Neural Networks (GNNs) to transform random initial points into low-discrepancy configurations, achieving state-of-the-an results by directly optimizing point coordinates (Rusch et al. 2024). Our work differs by framing the task as program synthesis problem rather than direct coordinate optimization.
related line of research focuses on optimizing the parameters of Sobol' sequences. The quality of Sobol' sequence is critically dependent set of initialization parameters known direction numbers. The direction numbers published by Joe & Kuo (2008) are a widely-used standard, de- rived from extensive computational search to find parameters that ensure high uniformity in two- dimensional projections by minimizing a quality measure known as the t-value. Subsequent work has focused further improving these parameters properties crucial for applications such as computer graphics (Bonneel et al. 2025). LatNetBuilder (eeucyer & Munger 2016) generalizes these ideas by searching the parameters of lattice rules and digital nets to minimize user-specified criteria, such as discrepancy bounds or projection-deeendent t-values.
guaranteeing quality for specific projection
Our approach is most closely related to the emerging paradigm of using Large Language Models (LLMS) for automated scientific discovery. Novikov et al. (2025) introduces AlphaEvolve, frame- work that combines LLMs with an evolutionary search, treating algorithm discovery as program evolution problem where the LLM functions as an intelligent mutation operator and receives feed- back from fitness function. This method has successfully discovered novel algorithms for funda- mental problems, from matrix multiplication (Fawzi et al. (0022) to open mathematical conjectures (Romera- Paredess et al. 2024). Our work is directly inspired by these principles, applying similar evolutionary loop to the specific mathematical challenges of discovering diw-dicreppancy sets and optimizing Sobol' sequences.
3

---

Published as a conference paper at ICLR 2026
4 METHODOLOGY
Our approach is based on the OpenEvolve framework, an open-source implementation of the princi- ples demonstrated by AlphaEvolve (Novikov et al. (025). frames the search for novel mathemat- ical constructs
evolutionary search over population of programs that generate them (Koza 1994). The LLM serves as a sophisticated mutation operator, guided by a fitness function.
The evolutionary loop Paredes et al. 2024) proceeds follows:
1. Initialization:
The process begins with population of "parent" programs. These pro- grams are code snippets in Python that generate a candidate solution. The initial population can be seeded with simple heuristics or well-known constructions.
2. Evaluation:
Each program in the population is executed, and its output is evaluated by fitness function. The fitness function returns a scalar score quantifying the quality of the solution (c.g., lower rQMC MSE or star discrepancy is better).
3. Selection and Prompting: High-performing programs are selected serve parents. A detailed prompt is then constructed for the LLM, including the parent program's is source code, its fitness score, and an instruction tasking the LLM with generating a variation that will improve upon the score. The prompt also includes code from other high- perforrmng "inspirations" to encourage crossover as well as guidance from the user (Appendix A).
4. Generation (Mutation): The LLM receives the prompt and generates a new, modified program. This is the core "mutation" step. The LLM's ability understand code syntax and semantics allows for complex and intelligent modifications.
5. Loop: The newly generated program is evaluated, its fitness is scored, and it is added to the population. The process then repeats, iteratively refining the population toward better solutions.
4.1 EXPERIMENTAL. SETUP
We used LLM-guided evolutionary search population of Python programs (population 60, four islands with occasional migration). Parents were chosen by fitness with occasional archive sampling; the LLM rewrote the mutation operator with temperature 0.7 and top-p 0.95. We did not train or fine-tune the LLM in any way; it is used as an off-the-shelf code rewriter. Each evolutionary run used roughly 2000 LLM calls and fixed budget for fitness evaluations. Full settings, including population structure, archive, migration, and program validation checks, are given Appendix A.
functions each generation. We used Google's Gemini 2.0 Flash as
4.2 DISCOVERY OF Low- NOW-DISCREPANCY POINT SETS
A two-phase strategy was employed to balance broad exploration with fine-tuned optimization.
Phase I: Direct Construction: The LLM was prompted to generate Python code that directly constructs an N N-point set in a d-dimensional space. The initial parent program in 2D implemented a simple shifted Fibonacci lattice (Appendix A, Listing 3 and in 3D implemented a scrambled Sobol' sequence (Appendix A, Listing This phase encouraged the LLM explore wide range of constructive heuristics.
Phase II: Iterative Optimization: After a sufficient number of iterations, the LLM prompted generate Python code that uses iterative optimization routines (e.g., s scipy. opt imiie..ii... ze. mim.mii ze) to refine an initial guess. This shifted the search from finding explicit constructions direct optimization the point coordinates.
Fitness Function: The fitness score was I+DE where Dyy is the star discrepancy of the generated point set.
We compute Dýy exactly by scanning extremal anchor boxes on the ceardinate-inducee grid with computational complexity (N44/+11). Implementation details in Appendix A.
4

---

Published as a conference paper at ICLR 2026
Initial Program
Direct Construction
Iterative Optimization
1.0
1.0
1.0
0.8
0.8
0.8
Yaxts
Yaais
0.4
0.4
0.4
0.2
0.2
0.2
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Xaxiss
Figure Visualization of N 16 point set generation two dimensions. (A) Initial shifted Fibonacci lattice (Discrepancy: 0.0962). (B) Best direct construction found in Phase (Discrepancy: 0.0924). (C) Final optimized point set from Phase II (Diserepancy: 0.0744), which is within 0.68% of the known optimal value of 0.0739.
Sobol
Rank-1 Lattice
Halton
MPMC
Hammersley
LLM
Star Discrepancy
Flbonacci
10-2
200
400
600 Number of Points (N)
800
1000
Figure 2: The Star Discrepancy Dy of Sobol', Halton, Hammersley, Fibonacci, Rank- Lantte-Lattice, MPMC (message passing Monte Carlo), and LLM-evolved sets for increasing number of points N 100.. 1020 in 2D.
4.3 DISCOVERY OF SOBOL' DIRECTION NUMBERS
Program Representation: The evolved programs Python functions that return list
of dictionaries. Each dictionary contains the Sobol' parameters (s, m_i) for a single dimension.
Initialization: The initial population contains the following implementation of the direc- tion numbers (Joe & Kuo 2008).
The primary fitness metric is 1+MSE where MSE is the mean squared error an rQMC estimate for a 32-dimensional Asian option price. The MSE calculated for N 8192 points and averaged over 1000 consistent randomizations (left matrix seram- ble followed by Cranley-Patterson random shift, i.e., LMS+shift) to ensure robustness and reproducibility. All rQMC comparisons are paired by using identical randomization seeds per method and N. The diffusion paths are constructed from [0, 1]d via standard time discretization geometric Brownian motion with equal timesteps; we report this mapping to clarify dffective-dimension effects.
Fitness Function:
5 EXPERIMENTS AND RESULTS
5

---

Published as a conference paper at ICLR 2026
Asian
Lookback
Barrier
Basket
Bermudan
Reduction MSE
-5
-10-
128
256
512
1024 2048 4096 8192
Number of Points (N)
Figure 3: The % reduction in MSE (rQMC integration over 10000 random scrambles and shifts) using Sobol' direction numbers found via LLM evolutionary search vs. those Joe & Kuo (2008). The reductions in MSE are averaged across all scenarios of that particular option (Appendix C).
Listing Directly Constructed 16 Point Set (N 16)
def construct_ste ()
2
A np. zeros (N, )
3
phi- (math. sqrt (5) /2
for i range (N)
5
A[1, 01= (1/math.mmaaeeeee sert (3))
6
A11, (]= 5+ (1)
7
return A
5.1 DISCOVERY OF LOW-DISCREPANCY POINT SETS
We first applied our two-phase strategy to the canonical problem of discovering N N-point sets in 2D and 3D with minimal star discrepancy. We illustrate the discovery process for 2D 16-point set (Fig. D. The initial program, a Fibonacci lattice, had a discrepancy 0.0962. After 243 iterations in the direct construction phase, new construction which consisted of fine-tuning optimal shifts to the Fibonacci lattice. After switching to the iterative optimization phase, the framework further refined the point set, achieving final discrepancy of 0.0744, which is within 0.68% of the known optimal value of 0.0739. The final program creates an initial guess, consisting of randomly jittered Fibonacci lattice, followed by sequential least squares programming (SLSQP) optimization loop with stochastic restarts (Listing 2).
found with discrepancy of 0.0924 (Listing [
We then benchmarked our method against Fibonacci. MPMC (Rusch al. 2024), and known opti- mal point sets (Clément et al. 22024) for N = 1... 100 [Table Our method successfully rediscov- ers the known optimal point sets for N < 10 and remains highly competitive for larger N. The most significant results came from searching for larger point sets where optimal solutions are not known. For 2D point sets with N 30, LLM evolutionary search discovered new configurations with lower star discrepancy than the best-known values from the literature. For instance, for N 100, our method found a point set with a discrepancy 0.01 50, substantial improvement over the previous best of 0.0188.. For N 100 we further compare the full evolutionary loop against single one-shot prompt and a multi-turn prompting baseline without evolution, as well as against a smaller Gemini Flash-Lite model. Over 16 seeds, evolutionary search with either LLM achieves consistently lower discrepancy with much smaller variance than the prompting-on baselines (Appendix B), showing that population- based evolution is more robust than repeated prompting.
We generated 2D point sets up to N = 1020 with lower star discrepancy than previously reported (Appendix B). In 3D, method matched the known optimal point sets for N 1,2,3, 6,

---

Published as a conference paper at ICLR 2026
N
Fibonacci
MPMC
LLM
Clément al.
1.0000
0.6180
0.6180
0.6180
0.6909
0.3660
0,3660
0.3660
0,5880
0.2847
0.2847
0.2847
0,4910
0.2500
0.2500
0.2500
0.3528
0.2000
0.2000
0.2000
0.3183
0.1692
0.1667
0.1667
0.2728
0.1508
0.1500
0.1500
0.2553
0.1354
0.1328
0.1328
0.2270
0.1240
0.1235
0.1235
10
0,2042
0.1124
0.1111
0.1111
0.1857
0.1058
0.1039
0.1030
0.1702
0.0975
0.0960
0.0952
13
0.1571
0.0908
0.0892
0.0889
14
0.1459
0.0853
0.0844
0.0837
15
0.1390
0.0794
0.0791
0.0782
0.1486
0.0768
0.0745
0.0739
17
0.1398
0.0731
0.0712
0.0699
18
0.1320
0.0699
0.0676
0.0666
19
0.1251
0.0668
0.0654
0.0634
0.1188
0,0640
0.0611
0.0604
30
0,0792
N/A
0,0438
0.0424
0.0638
0.0331
0.0332
50
0.0531
N/A
0.0278
0.0280
60
0.0442
0.0273
0.0234
0.0244
100
0.0275
0.0188
0.0150
0.0193
Table 2D Star Discrepancy Comparison for N 1 ... 100 between Fibonacci, MPMC (Message-
Passing Monte Carlo), LLM evolutionary search, and Clément et al. (provably optimal for N s 20). Best values are bolded.
Listing 2 Iteratively Optimized 2D Point Set (N 16)
def construct_steeee star ()
2
x zeros (N, ))
3 for i in range (N)
xii, 0]
np. rag.radoo.zand() /
x[i, 1] ((i 0.38196601125) 1)
rp. random, rand ()/ (2*N)
def discrepancy_wraperiss
7
points reshppp 2)
8
return star star_ discrepancy (points)
x0 x. x.
10
bounds = | (0.0, 1.0) (N 2)
11
best_result - None
12
best_disscrepancy
float (finf')
for range (25):
x0 20__rsttar x. flatten) ap. www.rando.normma.. (0, 0.01, N 2) 0.0, 1.0)
15
x0_restart =
16
result = minimize (disereanyyyy pepprppacyywwaappr, xt restatt,, method='
SLSOP' bounds=bounds, options=( maxiter'; 30000, ' 'foo'';
1c-15, 'iprint': 01)
17
discrepancy = discrepancy_
(result.x)
18
if discrepancy < best_di screpancy:
19
pest_discrepance
di screpancy
20
best_result result
potimized_poinne best_rresut..... .tt.trresstt..eesiliiee.. reshape (N, 2)
22
return optimized points
(Table 2) and provided explicit constructions that set new best-known durcriicrepancy benchmarks beyond the proven-optimal range for N (Appendix B).

---

Published as a conference paper at ICLR 2026
N
3
4
5
6
MPMC
0.6833
0.4239
0.3491
0.3071
0.2669
0.2371
0.2158
0.1993
LLM
0.6823
0.4239
0.3445
0,3042
0.2618
0.2326
0.2090
0.1937
Clément et al.
0.6823
0.4239
0.3445
0.3038
0.2618
0.2326
0.2090
0.1875
Table 2: 3D Star Diserepancy for N 1...8. Lower is better; best per N is bolded.
Training Example
Out-of-the-Money
At-the-Money
Sobol
LLM
p-value
Sobol
LLM
p-value
Sobol
LLM
p-value
32
0.248
0.252
0.976
0.192
0.198
0.976
0.325
0.335
0.976
64
0.0642
0.0665
0.998
0.0605
0.0631
0.998
0.0873
0,0908
0.998
128
0.0183
0.0192
1.00
0.0220
0.0231
1.00
0.0277
0.0282
1.00
256
0.00560
0.00610
1.00
0.00860
0.00890
0.00920
0.00980
1.00
512
0.00165
0.00161
0.284
0.00347
0.00331
0.0208
0.00325
0.00311
0.0208
1024
5,42c-04
5.24e-04
0,0548
0.00146
0.00131
3.5e-08
0.00123
0.00117
0.0304
2048
2.33c-04
2.25e-04
0.0199
6.1
5.35e-04
1.5e-14
5,50c-04
5.21e-04
0.0131
4096
9.88c-05
9.35c-05
0.0058
2.58c-04
2.43e-04
4.5e-04
2.40c-04
2.26e-04
0.0019
8192
4. 4.52c-05
4.10c-05
7.0c-06
1.17e-04
1.07e-04
1.0c-06
1.02e-04
9.52e-05
0.0100
In-the-Money
High Volatility
Low Volatility
Sobol
LLM
p-value
Sobol
LLM
p-value
Sobol
LLM
p-value
32
0.140
0.143
0.976
2.08
2.15
0.976
0.0238
0.0246
0.976
64
0.0367
0,0382
0.998
0.571
0.596
0.998
0.00660
0.00680
0.998
128
0.0106
0.0111
1.00
0.181
0.188
1.00
0.00217
0.00221
1.00
256
0.00308
0.00337
1.00
0.0598
0.0643
7.56e-04
8.03e-04
1.00
512
8.68c-04
8.59e-04
0.284
0.0205
0.0196
0.0073
2.85e-04
2.75e-04
0.0208
1024
2.61c-04
2.60e-04
0.258
0.00757
0.00706
3.4e-04
1.13c-04
1.09e-04
0.0717
2048
1.01c-04
1.00e-04
0.388
0.00315
0.00287
1.8e-66
5.21e-05
4.96e-05
0.0056
4096
4.04c-05
3.95c-05
0.218
0.00129
0.00123
0.0022
2.32e-05
2.20c-05
5.7c-04
8192
1.77c-05
1.68e-05
0.0047
5.66e-04
5.32c-04
0.0156
1.01e-05
9.30c-06
1.3e-04
Table 3: Mean Squared Error (MSE) and p-values for Asian option scenarios. The table compares the standard Sobol' sequence against sequence found via LLM evolutionary search (LLM). P- values are from one-sided Wilcoxon signed-rank test conducted over the 10000 randomizations with pairing by identical randomization seeds and are false discovery rate (FDR) corrected. P- values below 0.05 are bolded.
One natural concern is that Phase II might simply be rediscovering standard coordinate-wise local optimization, with the LLM contributing little beyond calling scipy. opt imi ze miniimzze. To test this, ran baselines that apply SLSQP directly Sobol' point sets, Fibonacci lattice, and best Phase I sets, all with the same objective and similar iteration budgets. These SLSQP baselines substantially improve over their respective seeds, but the LLM-evolved Phase programs still achieve lower star discrepancy in 22 of the 24 tested values of N with reductions of up to 15% relative to the best SLSQP-only baseline (Appendix B). This suggests the gains are not due to SLSQP alone, but to the LLM discovering nontrivial initializations and restarts that make SLSQP more effective.
5.2 IMPROVED DIRECTION NUMBERS
Having demonstrated the framework's 's ability to construct point sets, we next applied it to the discrete optimization problem of discovering improved Sobol' direction numbers. After several hundred evolutionary iterations, our L.LM evolutionary search routine discovered more performant set of parameters, focusing its modifications majority of the variance in Asian option pricing. Specifically, the parameters for dimensions 4, 5, and 6 were updated (Appendix C). All other dimensions (up 32) remained identical to the Joe & Kuo 2008 baseline.
the early dimensions, which are known to explain the vast
To validate these new direction numbers, they were benchmarked against the standard Joe & Kuo 22088) parameters across a suite of six Asian option pricing scenarios with varying parameters (Ap- pendix C). The true option price computed by taking the average of 1000 randomly scrambled
8

---

Published as a conference paper at ICLR 2026
Sobol sequences with N 2 points each. The rQMC MSE was evaluated over 10000 random seeds for = 32... 829 points. The direction numbers discovered by LLM evolutionary search produced a significantly lower integration MSE for larger sample sizes N z 512 under one-sided Wilcoxon signed-rank test and false discovery rate correction (Table
To disentangle the role of initialization from that of the LLM-guided search, we repeated the Sobol' experiment under three variants: (i) our main setup, warm-started from the Joe-Kuo parameters with Gemini Flash; (ii) the same pipeline with the smaller Gemini Flash-Lite model; and (iii) the same pipeline but initialized from random direction numbers rather than Joe-Kuo. Across six 32D Asian-style payoffs and sample sizes N, both Flash and Flash-Lite achieve similar MSE reductions over Joe- Kuo for N 2 2048, whereas the randomly initialized variant remains worse than Joe- Kuo for most N (Appendix C). This suggests that the evolutionary search is leveraging strong domain-specific prior given by Joe-Kuo, and refines it, rather than solving Sobol' design entirely from scratch.
We also compare our evolved direction numbers against Sobol' digital nets constructed by LATNET- BUILDER (L'ecuyer & Munger 2016), using its recommended random-CBC search and projection- dependent i-value criterion in d 32. For all six Asian-style payoffs and N S 1024, our direction numbers yield substantially lower rQMC MSE than the LatNetBuilder nets, often by factors of 2-10, with differences vanishing only for the largest N where all designs converge (Appendix C).
To ensure the evolved parameters were not merely overfitted to the Asian option's 's specific payoff structure, tested their generalizability on diverse suite of high-dimensional exotic options, in- cluding Lookback, Barrier, Basket, and Bermudan options (Appendices C, D). The evolved direction numbers demonstrated strong, generalizable performance, achieving significantly lower integration error across this wider range of financial instruments for larger sample sizes (N 2 512) with the sole exception of Barrier options. This suggests that the evolutionary search discovered a Sobol' sequence with fundamentally more robust and broadly applicable projection properties.
DISCUSSION
Using an LLM evolutionary framework, we generate 2D and 3D point sets with state-of-the-an star discrepancy and discover Sobol' direction numbers that lower rQMC integration error for high- dimensional exotic financial option pricing.
Recent state-of-the-art methods, such as the mathematical programming approach of Clément et al. (024) and the GNN-based MPMC framework of Rusch et al. (2024). are designed to construct point sets for fixed number of points N and dimensions d. In contrast, our approach generates the direction numbers that define a Sobol' sequence, offering several key advantages:
1. Extensibility to any N: Unlike fixed approaches where a new, computationally expensive optimization must be performed from scratch for each value of N, single, compact set of discovered direction numbers can be used to generate a high-quality point set for any desired number of points. This makes the solution immediately applicable to wide range of practical problems without requiring any re-computation.
2. Support for Progressive Integration: An integration can be performed with N points, and if more accuracy is needed, the next N points from the sequence can be added to refine the estimate while reusing all previous calculations. This fundamental capability that static point set generation methods inherently lack.
3. Easy Randomization: Unlike highly optimized, deterministic point sets, the LLM opti- mized Sobol' sequence can make use of standard randomization techniques, such Owen scrambling (Owen) 1995 1998), to obtain unbiased error estimates of rQMC.
4. High--Dimnnionas Applicability: The Sobol' framework is designed from the ground up for high-dimensional integration. By optimizing the parameters within this framework, our approach is directly applicable to problems such as the 32-dimensional option pricing benchmarks used in our tests, domain where direct coordinate optimization for a large N would be computationally intractable.
Other more recent systems, such as ShinkaEvolve (Lange et al. 2025). extend the same underlying AlphaEvolve loop with additional engineering features tailored to long-horizon code benchmarks,
9

---

Published as a conference paper at ICLR 2026
such richer parent sampling, novelty filtering, and advanced orchestration of program executions. In our setting, each candidate's 's fitness is a cheap, deterministic quantity (star discrepancy or rQMC MSE) computed by running short Python script, the core evolutionary mechanism provided by OpenEvolve is sufficient to explore the search space effectively. We deliberately avoided any RI.-style training or fine-tuning of the LLM: all improvements come from search in program space driven by an off-the-shelff model and the numerical fitness signal.
Taken logether with the broader suite of 32D lookback, basket, and Bermudan payoff' in Ap- pendix D. these results suggest that our evolved direction numbers are genuinely useful beyond the specific Asian option used in the fitness function: they systematically reduce MSE for range of smooth, path-dependent integrands. but do not help (and may slightly hurt) for highly discontinuous integrands such as barrier options. We view this as an instance of wroblem-dependen Sobol' design: the same pipeline could in principle be re-run with an objective that mixes several payoff types to trade off specialization and robustness.
Methodologically, our work does not introduce new evolutionary operators; our contribution is em- pirical and conceptual: we show that a generic LLM-guuded program search, applied almost out-of- the-box, is capable of rediscovering and subtly improving long-studied QMC constructions, and we analyze when these gains go beyond what can be achieved with standard local optimization alone. Empirically, our evolved Sobol' parameters work best for smooth, grands, and do not improve performance for highly discontinuous functions such as the payoff of barrier options.
variation inte-
One practical limitation of our approach is computational cost. For example, the full Sobol' search requires roughly 2000 LLM calls and on the order of x 1.6x 1020 32-dimensional QMC samples and payoff evaluations (Appendix A). On our CPU workstation this corresponds to about 96 wall-clock hours per Sobol' experiment. However, this cost is incurred once, offline: the resulting direction numbers can then be reused across many downstream Monte Carlo tasks.
We quantify significance using paired signed-rank tests and report all seeds; however, due com- pute limits we performed one evolutionary run per problem, which we treat as a limitation. Future work could explore alternative prompting techniques, optimization of evolutionary programming hyperparameters, and meta-model methods that generate optimal direction numbers for any given dimension point sets for any given dimension number of points N.
CONCLUSION
We have demonstrated the application of an LLM-driven evolutionary framework to tackle complex discovery problems in the field of low-discrepancy sets. Our method has discovered new 2D and 3D point sets with star discrepancy values lower than any previously published, setting new benchmarks in field of long-standing mathematical interest. Furthermore, it has produced novel Sobol' direc- tion numbers that improve the accuracy of rQMC integration for a variety of 32-dimensional finan- cial derivatives. This work strengthens the case for using LLMs core components in an automated scientific discovery process, capable of generating novel and valuable mathematical knowledge.
REFERENCES
Nicolas Bonneel, David Coeurjolly, Jean-Claude lehl, and Victor Ostromoukhov. Sobol' sequences with guaranteed-quality 2d projections. ACM Transactions on Graphics, 44, 2025. doi: 10. 3730821.
François Clément, Diederick Vermetten, Jacob de Nobel, Alexandre D. Jesus, Luís Paquete, and Carola Doerr. Computing star discrepancies with numerical black-box optimization algorithms. In Proceedings of the Genetic and Evolutionary Computation Conference, pp. 1330-1338. ACM, 2023. doi: 10. 11.1144/533313135904556
François Clément, Carola Doerr, Kathrin Klamroth, and Luis Paquete. Constructing optimal Loc star discrepancy sets. arXiv preprint
2024.
Josef Dick and Friedrich Pillichshammer. Digital Nets and Sequences. Cambridge University Press, 2010. doi: 10.1017/CB09778011777
10

---

Published as a conference paper at ICLR 2026
Alhussein Fawzi, Matej Balog, Aja Huang, Thomas Hubert, Bernardino Comera-Paredess Mo- hammadamin Barekatain, Alexander Novikov, Francisco J. R. Ruiz, Julian Schrittwieser, Grze- gorz Swirszcz, David Silver, Demis Hassabis, and Pushmeet Kohli. Discovering faster ma- trix multiplication algorithms with reinforcement learning. Nature, 610:47-53, 2022. doi:
10.10388441566--22-0511111
Gemini Team et al. Gemini: family of highly capable multimodal models. Report, 2025.
Google Technical
Paul Glasserman. Monte Carlo Methods in Financial Engineering, volume 53, Springer New York, 2003. doi: 10.1007/978-0-387-211111
Edmund Hlawka. Funktionen von beschränkter variation in der theorie der gleichverteilung. Annali di Matematica Pura ed Applicata, 54:325-333 333, 1961. doi: 10.1007/BF024155661
Stephen Joe and Frances Y. Kuo. Constructing Sobol sequences with better two-dimensional pro-
jections. SIAM Journal on Scientific Compuring. 30:2635 2654, 2008. doi: 10.1 10.1137/070709359,
F. Koksma. The theory of asymptotic distribution modulo one. Compositio Mathematica, 16: 1-22, 1964.
John R. Koza. Genetic programming as a means for programming computers by natural selection. Statistics and Computing, 4:87-112, 1994. doi: 10.1007/BF00115355.
Robert Tjarko Lange, Yuki Imajuku, and Edoardo Cetin. Shinkaevolve: Towards open-ended and sample-efficient program evolution. arXiv preprint arXiv:2509. /9349, 2025.
Pierre L'ecuyer and David Munger. Algorithm 958: Lattice builder: A general software tool for constructing rank-1 lattice rules. ACM Trans. Math. Softw. 42(2), May 2016. ISSN 0098-3500.
doi: 10.1 10.1145/2754929, RoLhttps: http://oo....gg9///// 1115//259114///74922223
Harald Niederreiter. Random number generation and quasi-Monte Carlo methods. Society for Industrial and Applied Mathematics, USA, 1992. ISBN 0898712955.
Alexander Novikov, Ngån Vü, Marvin Eisenberger, Emilien Dupont, Po-Sen Huang, Adam Zsolt Wagner, Sergey Shirobokov, Borislav Kozlovskii, Francisco J. R. Ruiz, Abbas Mehrabian, M. Pawan Kumar, Abigail See, Swarat Chaudhuri, George Holland, Alex Davies, Sebastian Nowozin, Pushmeet Kohli, and Matej Balog. AlphaEvolve: A coding agent for scientific and algorithmic discovery. arXiv preprint arXiv:2406. 13131, 2025.
Art B. Owen. Randomly permuted (t,s,))-nes and (t, s)-sequenees. Harald Niederreiter and Peter J.-S. Shiue (eds.), Monte Carlo and Quasi-Monte Carlo Methods in Scientific Comput- ing, volume 106 Lecture Notes in Statistics, pp. 299-317. Springer, 1995. doi: 10.1007/ 978-1-4612-25544 24.
Art B. Owen. Latin supercube sampling for very high- nigh-dimensional simulations. The Annals of Statistics, 26:619-640 1998. doi: 10.121 2/11////// (aaaoos///21111 0128144853.
Spassimir H. Paskov and Joseph F. Traub. Faster valuation of financial derivatives. The Journal of Portfolio Management, 22:1 13-123, 1995. doi: 10.3905/jpm. 1995. 409541.
Bernardino Romera-Paredes, Mohammadamin Barekatain, Alexander Novikov, Matej Balog, Pawan Kumar, Emilien Dupont, Francisco R. Ruiz, Jordan S. Ellenberg, Pengming Wang, Omar Fawzi, Pushmeet Kohli, and Alhussein Fawzi. Mathematical discoveries from program search with large language models. Nature, 625:468-477 475, 2024. doi: 10.1038/
$411866-033-09224-6
Konstantin Rusch, Nathan Kirk, Michael M. Bronstein, Christiane Lemieux, and Daniela Rus. Message-passing monte carlo: Generating low-discrepancy point sets via graph neural networks. the Proceedings the National Academy of Sciences, 121, 2024. doi: 10. 111003/paa..24099933.
Asankhaya Sharma. Openevolve: an open-source evolutionary coding agent, 2025. URL https 7/ wegithuue thub. http:///oooollloo///eeeeee/o/

---

Published as a conference paper at ICLR 2026
1. M. Sobol'. On the distribution of points in cube and the approximate evaluation of integrals. USSR Computational Mathematics and Mathematical Physics, 7:86-112, 1967. doi: 10.1016/ 0041-5553(677900444.9
12

---

Published as a conference paper at ICLR 2026
L.M Prompt for Direct Construction
You are an expert mathematician specializing in the construct on of QMC sampling points in a square. Your task to improve constructor function that directly outputs the position of 16 points on unit square ([0, 1] x [C, 11) a way that minimizes the star discrepancy. The star discrepancy is a measure of how uniformly is defined as the supremum of the absolute value of the difference between the fraction points and the area.
distributed the points are in the square.
Focus on designing position of each point (x, y) in the unit square, rather than iterative search algorithm.
explicit constructor that specifies the
It should output the position of each point (x, y) in the square [0, square.
10, 1l. 0.0 and 1.0 are included in the
Figure 4: The full prompt provided to the LLM for generating programs that directly construct point sets with minimum star discrepancy.
A IMPLEMENTATION AND COMPUTATION COMPLEXITY
A.1 INITIAL PROGRAMS
The initial parent program in 2D implemented simple shifted Fibonacci lattice (Listing) 3) and in 3D implemented scrambled Sobol' sequence (Listing
Listing 3 Initial Program for 2D Point Set Search
def construct_star ():
2
A-np. zeros (N, 2))
phi- (math .qrt (5) /2)
for
range (N):
5
/N
6
A AIi, 1 ((ipph 81+ (0.5/N) 81
7
return A
Listing Initial Program for 3D Point Set Search
1
def construct_star
2
A=Sobol (d=3, scramble=True, seed=42) random (n=N) return A
A.2 LLM PROMPTS
This appendix provides the full instructional text of the prompts provided to the Large Language Model (LLM) within the OpenEvolve framework. These LLM prompts were used to generate pro- grams to construct point sets directly (Fig. 4.. iteratively optimize point sets (Fig. 5). and directly
construct Sobol' direction numbers (Fig. 6
is important to note that the text shown in the figures below constitutes the instructional component of larger prompt. This text is dynamically combined with the source code of a "parent" program selected for mutation and often includes code from other nigh-peeforming "inspiration" " programs to encourage crossover of ideas.
13

---

Published as a conference paper at ICLR 2026
L.M Prompt for Iterative Optimization
You are an expert mathematician specializing in the const ruction of QMC sampling points in a 2D square. Your task is to improve a constructor function that finds the pesition of 16 points on a unit square (IO, 1] x 10, 11) in way that minimizes the star discrepancy,
The star discrepancy is a measure of how uniformly is defined as the supremum of the absolute value of the difference between the fraction of points and the area.
distributed the points are in the square.
Use scipy optimization routines such acipy. ooptiiize.minnnnne optimization routine and its initialization are critically important.
to fine-tune the construction. The
should output the position of each point (x, y) in the square [0, 11 10, 11. 0.0 and 1.0 are included in the square.
Figure 5: The full prompt provided the LLM for generating programs that iteratively optimize point sets to have minimum star discrepancy.
LLM Prompt for Sobol' Direction Numbers Search
expert mathematician specializing in the You ruction of QMC sampling points in square. Your task is to improve constructor function that directly outputs the direction numbers for dimensions 2 to 32 of a Sobol Sequence.
Your goal is to minimize the approximation error of 32 dimensional Asian option price. The dimensions 1, 2, and explain roughly 978 of the variance of the price.
The Sobol sequence is defined by a polynomial of degree s, with coefficients represented numbers ma for each dimension i. The direction numbers must be odd integers and within the specified range.
integer a, and direction
You must return list of 31 dictionaries for directions to 32, each containing the following keys:
"s" (int): The degree of the polynomial. 1 <= s <= 30 "a" (int): The coefficients of the polynomial, represented as an integer. 0 <= <
"m, (list The direction numbers for the Sobol sequence, represented as list of integers of length s. Each integer should be in the range [0, 2 and has to be odd.
Focus on designing an explicit constzuctor that specifies
these parameters, rather than iterative search algorithm.
The full prompt provided the LLM for generating program that directly specifies Sobol' direction numbers to have lower integration error for an Asian option.
Figure 6:
A.3 IMPLEMENTATION DETAILS
Our evolutionary algorithm treats each candidate solution as a full program (Python module) that im- plements either point-set generator (for the 2D star-discrepancy experiments) Sobol' direction-
14

---

Published as a conference paper at ICLR 2026
number constructor (for the QMC experiments). Each program exposes fixed entry point with a prescribed signature (e.g,, construct.poind (N) or construct_sobol .eeubbbeseeuunceee..... ( ) and is treated as black box: execute it in an isolated subprocess, validate its output, and then score it using the appropriate numerical objective.
Population structure and archive. We maintain population of P 60 candidate programs partitioned into four "islands" arranged in ring topology. Each island thus holds 15 programs. The search proceeds in asynchronous iterations: at each iteration we (i) select one parent program, (ii) mutate it via an LLM call to obtain single child, and (iii) evaluate and insert the child back into the perform migration:
island population. Every 25 iterations, sent clockwise and replaces the worst program at the destination. addition to the islands, we maintain global elite archive storing the top 25 programs
elite individual from each island
seen, ranked first by fitness and then lexicographically by (i) shorter program length and (ii) earlier discovery time. This archive used both for parent selection and for constructing the LLM context.
Parent selection and LLM-based mutation. At each iteration, we choose single parent for mutation. With probability 0.7, the parent is sampled from the elite archive; with probability 0.3, sampled from the union of all island populations. This biases search toward exploitation of high-performing programs while retaining exploration.
The mutation operator is call a Large Language Model (Google Gemini 2.0 Flash) configured as program rewriter. Rather than applying small, local edits, ask the LLM to produce complete rewrite of the core implementation functions, conditioned on several pieces of context:
the full source code of the parent program;
the source code of the three best-performing programs available ("inspiration programs"), which may come from any island or from the archive:
a system prompt specifying the required interfaces, validity constraints (e.g., dimensions, box [o, 1j odd direction numbers), and the target objective (minimizing 2D star discrep- ancy or randomized QMC MSE for the Asian option).
This setup acts as a soft form of multi-parent crossover: the LLM is free to borrow subroutines and structural ideas from multiple high-quality programs but must output single coherent child program. We use Gemini 2.0 Flash with temperature 0.7, top-p 0.95, and a maximum context length of 8192 tokens. In practice, most calls use around 4000 tokens (prompt plus completion). We do not tune these L.M hyperparameters per task or perform any RI. based finetuning or training on the LLM; all experiments share the same configuration.
Program execution and validity checks. Every newly generated program is executed in a dedi- cated subprocess via small wrapper script. The wrapper:
1. imports the program module;
2. calls the required entry function (e.g.,
());
3. serializes the returned object to a temporary file.
The main process then deserializes this object and performs a strict validity check. For Sobol' firection-number search, we require: 32 total dimensions; each dimension entry to contain fields s, a, and m; < < 30; 0 < a < 2 2 |m] = s; each direction number satisfies 1 < m, < 2" and is odd.
If any of these checks fails, if the program crashes, or if execution exceeds 600 second wall- clock timeout, the candidate is marked invalid and assigned zero fitness. For the 2D star discrepancy experiments, we analogously check that the program returns an N-point set in [0, 1]2 with the correct shape and bounds; wall-clock timeout of 600 seconds used.
Fitness evaluation. For valid programs, we compute fitness using the same numerical objectives in the main experiments:
15

---

Published as a conference paper at ICLR 2026
2D star diserepancy. Given an N-point set P [0, 1]2 returned by the program, com- pute the star discrepancy D* and define fitness as a monotone decreasing function of D (P): f I+D This preserves the ranking of candidate programs while providing scalar objective for evolution.
Sobol' direction numbers (randomized QMC MSE). For Sobol', the program outputs
list of 31 parameter triples (s, a, {m,}) defining dimensions ......32 32. To score candi- date A, we estimate the mean squared error (MSE) of randomized QMC estimator for the 32-dimensional Asian call option used throughout our QMC experiments. Each evaluation draws R 1000 independent randomizations (digital shifts and lower-triangllar scram- bling matrices) and, for each randomization, generates = 8192 Sobol' points in d' = 32 via compiled C++ generator. For each point set evaluate the Asian payoff and price estimator Č accumulate the squared errors (Ć Co)² relative the known ground-truth price Co. and define
1 f 1 EP (ŌW Co)2
Search budget and hardware. For each task (2D star discrepancy and Sobol' direction num- bers), at most 2000 evolutionary iterations, each evaluating exactly new candidate program. Including the initial population, this yields at most a few thousand distinct programs per run. In practice, nontrivial fraction of LLM-generated programs fail validation (e.g., due to syntax errors or invalid solutions) and are quickly discarded with zero fitness; the bulk of computation is spent on valid candidates that survive to be evaluated. All numerical evaluations (star discrepancy, Sobol' generation, QMC payoffs) are executed on CPU. Experiments were car- ried out on workstation with an AMD EPYC 7763 CPU and 64 GB of RAM. full Sobol' nuueerr run at this budget took approximately 96 hours of wall-clock time. We did not perform any hyperparameter tuning of the evolutionary algorithm; population size, archive size, migration schedule, parent selection probabilities, and LLM parameters are fixed across all ex- periments. Code and all generated point sets and Sobol' parameters are available at https: //qithu..com cot. wwt.:///w.ccccaaayyee/e/e//eeeea
A.4 COMPUTATIONAL COMPLEXITY
The total computational cost of our method has three main components: (i) LLM-based mutation; (ii) program execution and validation; and (iii) numerical fitness evaluation. We analyze each in turn.
LLM mutation. Let T prongi and T Tumpal denote the number of tokens in the input and output of one LLM call. Each mutation prompt consists of: a fixed system instruction block, the source code of the parent program, and the source of three top "inspiration" programs. If the typical program size is L lines of code, then both the prompt and the completion are O(L) in length, so that one mutation costs (Tpronpp + T ouppt) (L) tokens. With most CLLM 2000 LLM calls per run, the total LLM complexity per run is
(GI L).
In practice, we cap the context length at 8192 tokens, while observed calls average around 4000 tokens (prompt plus completion), the token budget scales linearly with the number of iterations.
Program execution and validation. Program execution cost per candidate scales with program length and the cost of constructing its output. Under design, programs generate their outputs symbolically (using explicit formulas and loops) rather than by performing expensive internal sam- pling. For the tasks considered:
in 2D star-discrepancy experiments, point generation costs O(N) per program;
in Sobol' numeee experiments, the program simply returns list of 31 parameter triples, which is (1).
Validity checks for Sobol' parameters are linear in the size of this list (constant with respect N), and validity checks for 2D point sets are linear in the number of points N. Compared to fitness evaluation, these costs are negligible.
16

---

Published as a conference paper at ICLR 2026
Numerical fitness evaluation. The dominant cost lies in computing the numerical objective for each valid program.
2D star discrepancy. For point set P with N points in = 2,3, we use a star discrep- ancy routine that scales exponentially with dimension ((N4/++1). Let (N) denote this cost. If at most K' distinct programs are evaluated in run, the total cost for star discrepancy is
OKK (as)
For ranges of N and d, this term is modest relative the Sobol' QMC experiments.
Sobol' direction numbers (QMC MSE). For Sobol', single fitness evaluation of didate requires = 1000 randomized QMC estimators, each with N = 8192 points and dimension = 32. The compiled C++ Sobol' generator produces N points in O(Nd) time per randomization, and the Asian payoff evaluation likewise costs ((Nd): each path requires cumulative sum d steps, exponentials, and payoff computation. Thus one fitness evaluation costs
O
With R = 1000, N = 8192, and d = 32, each evaluation processes R x N ~ 8.2 x 10t 32 92-dimensional points, corresponding to on the order of 2.6 10 scalar operations. If at most K < 2000 distinct Sobol' programs are evaluated, the total numerical cost per run is
O(K R-NN- M-N-dd d),
which in our setting corresponds to roughly K x Rx R 1.6 x 10tH 32-dimensional QMC samples and Asian payoffs per full evolutionary run.
Overall complexity. Combining the above, the total computational complexity per run can be summarized as
(CLLM L) +
O(K)
+ OKK
LLM mutations
program execution validation
fitness evaluations
where Ci = (NV) for 2D discrepancy and Cn: = RNd for Sobol'. all our Sobol' experiments the K KRN term dominates, while for 2D discrepancy the LLM term and the discrepancy term of comparable size. Crucially, this entire cost is incurred once, offline. Using the resulting point sets Sobol' parameters in downstream applications has the same asymptotic cost as using a standard hand-designed formula: generating N Sobol' points is still ((Nd), and evaluating QMC estimators with them is unchanged.
B GENERATED POINT SETS
This appendix provides the full numerical results for the star discrepancy (Dy) values of the point sets discovered by LLM evolutionary search. These tables serve comprehensive record of the performance of our method, establishing new state-of-the-art benchmarks and offering detailed comparison against established fow-discrepancy construction methods.
B.1 STAR DISCREPANCY EVALUATION
We evaluate Dj by enumerating all anchored boxes whose upper corners lie on the Cartesian prod- of per-dimension grids G {unique ziz} [11, with the lower corner fixed at the origin. For each candidate corner y G1 X ... X Gd, we compute its volume II, Vj and two point counts: the number of samples dominated by (component-wise <) and the same count with points on any upper face removed. The local discrepancy is the maximum of the two absolute differences from the volume. capturing the half-open convention on anchored boxes. Inputs are clipped to [0, 1]d, grids are built from unique coordinates (plus the terminal 1.0).
B.2 NEW BENCHMARKS IN THREE DIMENSIONS
We report high-quality constructions for N from 16 in three dimensions (Table 4
17

---

Published as a conference paper at ICLR 2026
N
LLM Discrepancy (Dt)
0.1758
10
0.1652
11
0,1551
12
0,1483
13
0.1402
14
0.1337
0.1275
16
0.1207
Table 4: The Star Diserepancy Values for N
3D point sets found via LLM evolutionary search.
N=30
N=50
N=80
0.06
N=40
N=60
N= 100
0.05
Star Discrepancy
0.04
0.03
0.02
0.01
250
500
750
1000
1250
1500
1750
2000
Iteration
Figure 7: Best-so-far star discrepancy 2D for E {30, 40, 50, 60, 80, 100} as a function of iteration. Each curve corresponds to single evolutionary run; the value iteration t is the minimum discrepancy observed up to that point.
B.3 OPTIMIZATION TRAJECTORIES IN 2D
To better understand how the evolutionary loop improves the point sets over time, we record the best-so-far star discrepancy after each iteration for several values of N 2D. Figure 7 shows the resulting trajectories for N e {30, 40, 50, 60, 80, 100} over the full 2000-iteration budget.
Across all N, we observe a characteristic two-phase behavior. First, the discrepancy improves rapidly during the early iterations, roughly the first 200-400 steps, dropping sharply from the initial constructions.
After this initial phase, the curves enter slower regime of incremental refinement. They continue to decrease, but with smaller step sizes and less frequent improvements. Occasional late "jumps"- example, for N 80 around iteration 1 1100 and for N 30 around iteration 1200 that useful innovations can still be discovered well into the run. However, these gains are modest compared to those obtained early on.
B.4 EXTENDED PERFORMANCE IN TWO DIMENSIONS
In two dimensions, while the problem is better understood than in 3D, optimal solutions for larger point sets (N > 21) are not known. We provide detailed comparison of the 2D star discrepancy values achieved by LLM evolutionary optimization against several key baselines for ranging from 140 1020 (Table 5
18

---

Published as a conference paper at ICLR 2026
N
Halton
Sobol'
Hammersley
Fibonacci
MPMC
LLM
140
0.03686
0.02794
0.02701
0.01870
0.01373
0.01151
180
0.03200
0.02300
0.02283
0.01454
0.01147
0.01058
220
0.02323
0.01924
0.01868
0.01190
0.00975
0.00891
260
0.02062
0.01546
0.01581
0.01063
0.00843
0,00831
300
0.01994
0.01341
0.01424
0.00972
0.00752
0.00741
340
0.01950
0.01353
0.01266
0.00858
0.00710
0.00696
380
0.01659
0.01167
0.01170
0.00768
0.00695
0.00638
420
0.01617
0.01194
0.01058
0.00694
0.00584
0.00563
460
0.01279
0.00972
0.00966
0.00634
0.00540
0,00471
500
0,01172
0.00848
0,00889
0.00583
0.00518
0.00509
540
0,01101
0.00967
0.00823
0.00540
0.00488
0.00407
580
0.01261
0.00956
0.00766
0.00503
0.00476
0.00377
620
0.01140
0.00782
0.00744
0.00470
0.00454
0.00386
660
0.01074
0.00956
0.00699
0.00463
0.00437
0.00360
700
0.00933
0.00865
0.00682
0.00437
0.00416
0.00379
740
0.00891
0.00709
0.00646
0.00435
0.00397
0.00346
780
0.00870
0.00662
0.00612
0.00412
0.00376
0.00310
820
0.00881
0.00678
0.00583
0.00392
0,00360
0.00294
860
0.00914
0.00604
0,00556
0.00374
0.00356
0.00316
900
0.00751
0.00561
0.00531
0.00357
0.00319
0.00296
940
0.00785
0.00481
0.00508
0.00342
0.00321
0.00289
980
0.00768
0.00558
0.00487
0.00328
0,00308
0.00280
1020
0.00777
0.00395
0.00468
0.00315
0.00303
0.00245
Table 5: 2D Star Discrepancy Comparison for N 140. This table compares the performance of LLM evolutionary search against classical low-discrepancy sequences and the state-of-the-an MPMC method for larger point sets where optimal solutions are not known. Lower values indicate better uniformity. The best result in each row is bolded.
Consistent Outperformanee: Across every single value of N tested, the point sets dis- covered by LLM evolutionary search achieve a lower star discrepancy than all other meth- ods, including classical sequences (Halton, Sobol', Hammersley, Fibonacci) and the recent state-of-the-art MPMC method.
Significant Improvement Margin: The performance gap not trivial. For example, = 140, the LLM discrepancy of 0.01151 is approximately 16% lower than the next-best method (MPMC at 0.01373) and over 57% lower than the widely used Halton sequence. At N 1020, the LLM discrepancy of 0.00245 is nearly 20% better than MPMC's 0.00303.
The results demonstrate that the LLM evolutionary search's 's ability to find superior configurations is not limited to specific range of N, but holds consistently as the number of points increases. This suggests that the two-phase discovery strategy (direct construction followed by iterative optimization) is effective at navigating the increasingly complex search space associated with larger point sets.
Superior Scaling:
B.5 PIASEEIII LOCAL OPTIMIZATION
To determine whether Phase II improvements could be attributable solely the use of standard local optimizer, we constructed a suite of baselines in 2D that used purely local optimization, specif- ically sequential least squares programming (SLSQP) which the optimizer predominantly used by the LLM constructed programs. For each N E {2,... 20, 30, 40, 50, 60, 100} we consider three initial point sets: the Sobol' sequence, a simple Fibonacei construction, and the best point set re- turned by the LLM's Phase constructive search. From each of these initializations, we then run single SLSQP local optimization directly on the point coordinates, using the same star-discrepancy objective as in the main experiments and box constraints [0, 1|2 for all points. The resulting sets are denoted SLSQP, Fibonacci+SLSQN and Phase ++SLSQP. Importantly, these baselines use the same objective and comparable iteration budget to Phase II, but without any -es-indgg program structure (no LLM-chosen restarts, parameterizations, or heuristics). For comparison, we also report the M.L--vevlled point sets from Phase II ("LIM"').
19

---

Published as a conference paper at ICLR 2026
Sobol
Fibonacci SLSQP
Phase I
N
Sobol
SLSQP
Fibonacci
Phase
SLSQP
LLM
2
0.7500
0.3660
0.5365
0.3660
0.4375
0.3660
0.3660
0.6250
0.3333
0,4850
0.3333
0.3375
0.3333
0.2847
0.4375
0.2500
0.3637
0.2500
0.2856
0.2500
0.2500
5
0,4375
0.2000
0.2910
0.2000
0,2500
0,2500
0.2000
0.2917
0.1833
0.2836
0.1795
0.2078
0.1771
0.1667
7
0.3393
0.1571
0.2431
0.1587
0.1870
0.1518
0.1500
0.3125
0.1495
0.2127
0.1389
0.1697
0.1429
0.1328
0.3264
0,1318
0.1891
0.1270
0,1551
0.1270
0.1235
10
0,2906
0.1318
0.1702
0,1143
0.1375
0.1172
0,1111
11
0.2088
0.1099
0.1547
0.1074
0.1255
0.1095
0.1999
12
0.2240
0.1303
0.1418
0.0985
0.1258
0.0979
0.0960
0.2404 13
0.1176
0.1309
0.0909
0.1161
0.0918
0.0892
14
0.1970
0,1092
0.1295
0.0850
0.1075
0.0857
0.0844
0.2094 15
0.1020
0.1209
0.0793
0.1012
0.0804
0.0791
16
0.1719
0.0945
0.1264
0.0743
0.0962
0.0757
0.0745
17
0.1976
0.0937
0.1190
0.0732
0.0903
0.0735
0.0712
18
0.1649
0.0980
0.1124
0.0691
0.0885
0.0676
0.0676
19
0.1357
0.0969
0.1065
0.0673
0.0850
0.0665
0.0654
20
0.1313
0.0824
0.1011
0,0642
0.0818
0.0633
0.0611
30
0.0690
0.0613
0.0674
0.0442
0.0595
0.0440
0.0438
40
0,0836
0.0481
0.0579
0,0356
0.0457
0.0365
0.0331
50
0.0745
0.0379
0.0463
0,0281
0.0355
0.0278
0.0278
60
0.0484
0.0349
0.0386
0.0236
0.0320
0.0253
0.0234
100
0.0423
0.0227
0.0245
0.0153
0.0212
0.0161
0.0150
Table 6: 2D star discrepancy comparison for N E (2,..., 20, 30, 40, 50, 60, 100) between Sobol, Fibonacci, Phase (with and without SLSQP refinement), and LLM search. Best values are bolded.
reppots the resulting star discrepancies. As expected, adding SLSQP on top of classical con- structions already yields large improvement over the Sobol, Fibonacci, and Phase sets, confirm- ing that local optimization is a strong baseline for this problem. However, the LLM-evolved Phase II sets consistently achieve even lower star discrepancies. In 22 out of the 24 values of N in Table the LLM achieved strictly lower star discrepancy than every SLSQP baseline (Sobol+SLSQQ Fi- bonacci+$LSQPP, Phase 1+SLSQP). The largest relative improvements occur at small and moderate N (e.g., 0.3333 > 0.2847 at N = 3, reduction of about 15%. and 0.0356 - 0.0331 at N = 40, roughly 7%). In remaining two cases (N 16 and N = 50), the best SciPy baseline is better by less than 2 x 10-4 in absolute discrepancy.
These results show that while local optimization is essential for obtaining strong point sets, the L.M-- guided Phase II search still provides measurable advantage over just local optimization starting from strong human-designed seeds.
B.6 PROMPTING VS EVOLUTIONARY OPTIMIZATION AND LLM SIZE
Method
Mean
Std. dev.
Min
Max
singleturn
0,02117
0.00012
0.02090
0.02135
multiturn
0.01852
0.00254
0.01479
0.02098
L.LM
0.01541
0.00065
0.01492
0.01756
M.Ie-lite
0.01540
0.00064
0.01492
0.01690
Table 7: Ablation at N 100 (2D): star discrepancy over 16 seeds for single-shot prompting (singleturn), multi-turn prompting without evolution multiturn)) and evolutionary search with the Gemini Flash model (LLM) or the smaller Flash-Lite model (LLM-lite)
We next isolate the contributions of (i) the evolutionary search loop versus plain iterative prompting, and (ii) the choice of L.M.M At N 100 in 2D, we ran 16 independent seeds for four variants:
20

---

Published as a conference paper at ICLR 2026
Performance for N=100
0.021
Star Discrepancy
0.020
0.018
0.017
0.016
0ge0 Evolutionary (Flash-Lite)
0.015
Single-turn
Multi-turn
Evolutionary (Flash) Optimization Protocol
Figure 8: Per-seed star discrepancy N 100 in 2D for single-turn and multi-turn prompting, and evolutionary search with Gemini Flash and Flash-Lite. Each point one run with a different seed.
singleturn: a single one-shot call to the main LLM (no iteration, no evolution);
multiturn: repeated prompting of the LLM for up 2000 calls, but without evolutionary operators (no population, no selection, no mutation); this is an iterative prompting baseline;
LLM: the full evolutionary pipeline with the main Gemini Flash model, including population-based evolutionary search (selection, mutation, and archiving);
LLM-lite: the same evolutionary pipeline as LLM, but using a smaller Flash-Lite model as the LLM.
All variants optimize the same 2D star-discrepe dicredisrrppancy objective as in the main experiments. Table summarizes the distribution of the best star discrepancy for N = 100 points over the 16 seeds for each method, and Figure visualizes the per-seed outcomes.
First, the single-shot baseline insufficient to reach the low-discrepancy regime we target. Sec- ond, simply calling the LLM multiple times without an evolutionary mechanism (multi titurn) substantially improves over single-turn prompting (mean 0.0185 vs. 0.0212) and occasionally dis- covers very strong solutions (best seed 0.01479). However, this baseline exhibits high variability across seeds (standard deviation Rs 2.5 x 10 10-3) 3), with some runs remaining close to the single-turn performance. This variability is clearly visible in Figure
Third, the full evolutionary search variants (LLM and LIM-lite) achieve consistently lower dis- crepancies average (=s 0.0154) with much smaller variance (Rr 6.5 10-4). Their best (0.01492) are within 1.3 x 10-4 of the best mu titurn run, but the evolutionary runs are far more robust: almost all seeds converge to strong solutions, seen in the tight clusters on the right of Figure 8
The comparison between LLM and LLM LLM-lite further indicates that the overall evolution- ary search pipeline remains effective even when replacing the main LLM with smaller Flash-Lite model; their means and variances are nearly identical.
C GENERATED SOBOL' DIRECTION NUMBERS
S
A
M;
4
3
135
5
3
2
137
6
4
1
1137
Table Sobol' direction number parameters updated by LLMM evolutionary search. the dimen- sion, S is the polynomial degree, A is the polynomial's 's coefficients, and M, are the initial direction numbers. All other dimensions remained unchanged.
21

---

Published as a conference paper at ICLR 2026
We provide the direction numbers discovered by the LLM evolutionary search routine. Only the parameters for dimensions 4, 5, and 6 were updated (Table 8
C.I COMPARISON WITH LATNETBUILDER DIGITAL NETS
We compare Sobol' direction numbers found via evolutionary LM search against digital nets con- structed with LATNETBUIL a state-of-the-an QMC design tool. We use the recommended random-CBC search (2000 iterations) in LatNetBuilder with a projection-dependent t-value based criterion in = 32 dimensions for 2 points:
Latretbuilder -t net -c sobol -s N -d 32 -e MB:dom-CBC:2000 -f projdep:t value -q inf -W order- depeer--dpendett:::::0
For each resulting Sobol digital net and each number of points N E {32, 64..... 8192), we con- struct randomized QMC estimators for the six 32-dimensional Asian-style aption-pricing payoffs considered in the training example, sat-of-the-moncy at-the-money. nn-the-money, high volatility, and low volatility. Table reports, for each payoff' and N, the empirical MSE for LatNetBuilder vs. our implementation (LL.M), together with a p-value from a one-sided Wilcoxon signed-rank test with false discovery rate correction across payoff's for each N; significant p-values (p 0.05) bolded.
For small to moderate sample sizes (N < 1024), the LLM discovered Sobol' sequence consistently and substantially outperforms the LatNetBuilder Sobol' digital nets, often reducing MSE by factors between 2 and 10x, with extremely small p-values across all scenarios. As N increases further, the MSEs of the two constructions converge. In this regime differences are small in absolute terms and typically not statistically significant. Overall, the Sobol' direction numbers found by our imple- mentation are competitive with, and usually improves upon, Sobol' direction numbers optimized by a dedicated QMC design tool over a wide range of sample sizes.
C.2 ROLE OF LLM SIZE AND INITIALIZATION FOR SOBOL' DIRECTION NUMBERS
To better understand the role of LLM size and the importance of initialization in our direction- number search, we compare three LLM variants against the standard Joe-Kuo Sobol' sequence on the six 32D Asian-style options used in our main randomized QMC experiments. For each variant, we run 16 independent evolutionary searches with different random seeds and report the mean + standard deviation of the resulting randomized QMC MSE. The direction numbers used for evaluation are exactly those optimized on the training example; they are not retuned for the other payoffs. The results are summarized in Table TO
The first two variants differ only in the underlying LLM: LLM uses Gemini Flash, while LLM- lite replaces it with Gemini's 's smaller Flash-Lite model. Across all scenarios and sample sizes N, their MSEs are very similar, with differences typically well within one standard deviation. On the training example, for instance, both variants have slightly worse MSE than Joe-Kuo for small N but achieve lower MSE for larger N 2 2048. A similar pattern holds for the out-of-the-monee at-the- money, and volatilityshhitted payoffs: for N 2 2,000, both LLM-based variants match modestly improve on the Joe-Kuo baseline across almost all scenarios, despite being optimized only on the training payoff. This suggests that the overall evolutionary pipeline is not overly sensitive LLM size; a smaller model can still discover
direction numbers. Furthermore, both LLM and LLM-lite have relatively low standard deviations (typically less than 25% of the mean for large N), reflecting their stability across evolutionary runs.
The third variant, LLM (random init), replaces the Joe-Kuo direction numbers with random ini- tialization in the same parameter space, keeping the rest of the OpenEvolve pipeline unchanged. This variant is consistently worse than standard Sobol' at small and moderate (often by factors of 2x-5x in MSE), but approaches baseline performance at the largest N (e.g., N 8192). In several scenarios, such as at-the-money and high volatility, it remains worse than Kuo- These results indicate that, under realistic compute budget, evolutionary search alone is not sufficient to recover high-quality Sobol' parameters from scratch; a strong initialization such as Joc-Kuo is a crucial inductive bias.
22

---

Published as a conference paper at ICLR 2026
D OPTION SCENARIOS
To ensure a rigorous evaluation of the Sobol' direction numbers discovered, we designed a com- prehensive suite of benchmark scenarios. This suite includes the primary optimization target (the Asian option, see Table [1] as well as a diverse set of exotic options known to be challenging for Quasi-Monte Carlo integration. The purpose of this suite twofold: first, confirm superior per- formance on the target problem class, and second, to test for generalizability and ensure that the evolved parameters were not merely overfitted the specific payoff structure of the Asian option.
The configuration parameters for all tested scenarios are consolidated in Table 12 For all options, the time to expiration (T) was set 1.0 year and the risk-free interest rate (r) was 0.05. The under- lying asset prices are assumed to follow geometric Brownian motion.
D.1 ASIAN OPTIONS
An Asian option a path-dependent exotic option whose payoff is determined by the average price of the underlying asset option, which only depends on the asset price at expiration. The averaging feature reduces volatility and makes the option generally cheaper than its European counterpart. Because its value depends on the entire price path, pricing it requires simulating all 32 time steps, making it an excellent candidate for QMC methods.
pre-set period of time. This is in contrast to standard European
D.2 LOOKBACK OPTIONS
A lookback option is another path-dependent option whose payoff is determined by the maximum or minimum price the underlying asset over the option's life. The scenarios tested here are for floating strike lookback call option, whose payoff at expiration is the difference between the final asset price and the minimum price achieved (ST - Smsn).
D.3 BARRIER OPTIONS
barrier option is a path-depended option that is either activated ("knocks-in") extinguished o"knocks-out") if the underlying asset price crosses a predetermined "barrier" level. This feature introduces significant discontinuity in the payoff function.
D.4 BASKET OPTIONS
A basket option's payoff depends the value of portfolio or "basket" of multiple underlying assets. It is an inherently high-dimensional problem where the correlation (p) between the assets is critical parameter. For these tests, assume uniform initial price of So 100.00 for all 32 assets in the basket.
D.5 BERMUDAN OPTIONS
Bermudan option is a hybrid between European option (exercisable only expiration) and an American option (exercisable at any time). It can be exercised on discrete set of pre-specified dates. Pricing a Bermudan option is a highly complex problem that requires solving dynamic programming problem determine the optimal exercise strategy.
E RQMC INTEGRATION RESULTS
This appendix provides the detailed results for the primary benchmark (Asian option) and the gen- cralizability tests (exotic options). Each table compares the performance of the standard Sobol' sequence (Joe & Kuo) against the direction numbers discovered by LLM evolutionary search. We report the Mean Squared Error (MSE) and its constituent parts, Squared Bias and Variance, for an increasing number of points (N). We consistently found lower variance and MSE for all option types for all 2 512 with two exceptions: close barrier option (N = 2048) and high correlation basket option (N 2048).
23

---

Published as a conference paper at ICLR 2026
Scenario
N
LatNet MSE
LLM MSE
p-value
Training Example
32
3.472e-01
2.523c-01
1.97e-44
64
1.308e-01
6.649c-02
1.55e-177
128
7.057e-02
1.922c-02
0.00c+00
256
1889c-02
6.108e-03
0.00e+00
512
1.384c-02
1.614e003
0.00e+00
1024
9.989e-03
5.244e-04
0.00e+00
2048
2.753e-04
2.254c-04
4.15e-15
4096
1.213c-04
9.346c-05
9.37e-27
8192
4.286c-05
4.104e-05
1.63c-01
Out-of-the-Monee Money
32
2.532e-01
1.979-001
7.18e-25
64
1.002-001
6.309e-02
6.55e-77
128
5.247c-02
2.313e-02
8.85e-261
256
1.597c-02
8.893e-03
1.46e-148
512
1.087e-02
3.311c-03
0.00e+00
1024
6.971e-03
1.313e-03
0.00c+00
2048
5.971e-04
5.353c-04
3,29c-07
4096
2.459c-04
2.432e004
1.77c-01
8192
8.641e-05
1.066c-04
1.00c+00
At-the-Money
32
4.503e-01
3.346c-01
1.46e-42
64
1.691c-01
9.079e-02
1.74e-152
128
9.224c-02
2.822e-02
0.00e+00
256
2.413e-02
9.815e-03
0.00e+00
512
1.794c-02
3.106e-03
0.00c+00
1024
1.269c-02
1.168e-03
0.00c+00
2048
5.418c-04
5.212e-04
3.04c-01
4096
2.355c-04
2.261e-04
1.18c-01
8192
9.41le-05
9.523e-05
1.00e+00
In-the-Money
32
1.925c-01
1.428-001
1.95e-37
64
7,452c-02
3.818e-02
2.27e-168
128
3.849e-02
1.110c-02
0.00e+00
256
1.037e-02
3.373c-03
0.00c+00
512
7.949e-03
8.593e-04
0.00c+00
1024
5,697c-03
2.601e-04
0.00e+00
2048
1.315c-04
1.000e-04
2.91e-33
4096
6.003e-05
3.949e-05
2.16e-72
8192
1.693e-05
1.683e-05
6.02e-01
High Volatility
32
2,902e+00
2.149e+00
1.25e-40
64
1.13e+000
5.958e-01
1.74e-152
128
6.053c-01
1.877e-01
0.00c+00
256
1.647c-01
6.425e-02
0.00e+00
512
1.235c-01
1.959e-02
0.00e+00
1024
8.672c-02
7.060e-03
0.00e+00
2048
3.225e-03
2.867c-03
1.75e-04
4096
1.365e-03
1.233e-03
4.04c-05
8192
4.255e-04
5.321e-04
1.00++00
Low Volatility
32
3.253e-02
2.455c022
2.93e-38
64
1.200e-02
6.833c-03
2.93e-129
6.584e-03
2.205c-03
0.00c+00
256
1.726c-03
8.026e-04
6.07e-252
512
1.251c-03
2.747e-04
0.00e+00
1024
8.675c-04
1.088e-04
0.00e+00
2048
5.024c-05
4.961e-05
5.67c-01
4096
2.207e-05
2.196e-05
6.46i-011
8192
9.317c-06
9.300e-06
1.00e+00
Table 9: Mean-squared error (MSE) of randomized QMC estimators for six 32-dimensional Asian- style option payoffs, comparing digital nets constructed with LATNETBUILDE to those optimized by LLM evolutionary search, for N E (32, 64,..., 8192). Lowest MSE in each row bolded. The last column reports FDR-corrected one-sided Wilcoxon p-values; significant differences (p < 0.05) bolded.
24

---

Published as a conference paper at ICLR 2026
LLM MSE
Scenario
N
Sobol MSE
LLM MSE
LLM-lite MSE
Random Init
Training
2.48c-01
3.12e-01 -
6.20c-01 - 1.81e-01
Example
6.42e-02
1.05--01 8.77e-02
1.11e-01 5.53e-022
2.85e-01 + 1.04e-01
128
1.83e-02
3.39c-02 + 3.80e-02
3.75e-02 2.233e--0
1.34c-01 5.56e-02
256
5.55c-03
1.26e-02 + 1.85e-02
+ 8.09e-03
5.90e-02 2.24e-02
512
1.65-033
4.22e-03
= 6.52e-03
2.64e-03 = --63e-03
1.95e-02 9.69e-03
1024
5.42e-04
1.66--33 3.39e-03
7.97e-04
0.39
7.28e-03 4.31e-03
2048
2,33e-04
2.26e-04:11 3.91e-66
2.43e-04 + 1.65e-05
2.52e-03
2.75e-00
4096
9.88c-05
9.72c-05 + 3.02c-06
1.01c-M4 7.700-06
3.49c-04+4..88e 4.03e-00
8192
4.52c-05
4.18c-05 1.14c-06
4.25e-05
1.65e-06
4.53e-05 2.57e-06
Out-of-the
1.92-001
2.28e-01 8.22e-02
2.32e-01
1:: 4.92e-00
3.95e-01
1.0le-01
Money
64
6.05e-02
8.39e-02
1 4.89e-02
8.83e-02 3.15e-00
1.83e-01
1 5.91e-02
128
2.20c-02
3.06c-02
= 2.07c-02
3.33c-02 1.42cc002
+..00--33300000
256
8.60c-03
1.25c-02 = 1.05c-02
4.40c-03
+.75e-22+1.222-0
3,47c-03
4.93c-03 3.82c-03
+27c-02+4.999)
1.46--03
2.02e-03 + 1.94e-03
1.62e-03
3.36e-04
4.83e-03 2.24--03
2048
6.19e-04
5.44c-04 + 1.42e-05
5.93e-04
::4044:444 4.38e-005
1.74e-03 1.43e-03
4096
2.58c-04
2.48e-04 + 6.55e-06
2.57e-04 2.43e-005
3.59e-04
1.84e-01 1.84e-04
8192
1.17e-04
1.12e-04 = 4.3
1.10-04 = 1.18e-05
1.10e-04 1.07e-05
At-the
3.25c-01
4.04c-01 + 1.82e-01
4.12c-01 1.12e-01
7.76e-01
2.20e-01
Money
64
8.73e-02
1.38--01 + 1.08e-01
1.47e-01 + 6.94e-002
3.56e-01
+ 11227-001
128
2.77e-02
4.56e-02
+ 025
1.64e-01
+ 6.51e-02
256
9,23e-03
1.77e-02 2.26e-02
1577-02
7.20e-02 2.65e-02
3,25e-03
6.36e-03 8.09e-03
2.33e-02 1 1.10e-02
1024
1,23c-03
2.54c-03 4.09c-03
8.69c-03 4.87c-03
2048
5.50c-04
5.14c-04 + 1.52c-05
+.36--0422.66600
3.00c-03 + 3.08c-03
4096
2.40c-04
2.30c-04 7,48--06
8192
1.02-004
9.95e-05 3.97e-06
6.50e-05-1.188--0
1.03e-04 7.57e-06
In-the Money
I. 140c-01
1776-001
3.45c-01 1.02c-01
64
3.67e-02
5.85c-02 4.74c-02
1.61c-01 + 6.00c-02
128
1.06c-02
1.93c-02 = 2.16c-02
1.
4.33-0
256
3.08e-03
6.98e-03 + 1.03e-02
7.14e-03 4.700e-0
3.47e-02 #: 1.33e-02
8.68e-04
2.28e-03 + 3.56e-03
1.06-03
1.45e-03
1.06e-03
1.21e-02 + 6.13e-03
1024
2.61e-04
8.97c-04 + 1.94e-03
4..4-0444.24e-0
4.24e-004
4.61e-03 + 2.777-00
2048
1.0le-04
1.01e-04
= 1.70e-06
1.10-04411.022 1.02e-05
1.64e-03 + 11777e-03
4096
4.04e-05
4.18e-05 1.35e-06
4.36e-05 4.54e-06
2.25e-04 2.93e-04
8192
1.77e-05
1.71e-05 4.83e-07
1.72--05-8..41e
1.84-05 1.43e-06
High
32
2.08e+00
2.60e+00
1.19e00 1.19e+00
264e+00
3.0e=400
Volatility
64
5.71e-01
9.00e-01 + 7.05e-01
2.38e+00 +: 8.77e-01
128
1881-01
3.05e-01
+:2 3.09e-01
3.36e-01
2.07e-01
1.14+000 + 4.65e-01
256
5.98c-02
1.17c-01
1.51c-01
6.660-066666600
5.02c-01 = 1.86c-01
512
2.05c-02
4.14c-02 5.37c-02
1.70c-01 = 8.15c-02
1024
7.56c-03
1.66c-02 2.81c-02
9.77c-03=
5..1lc-03
6,43c-02 + 3.65e-02
2048
3.15e-03
2.87e-03 7.08e-05
114e-003:12...e
2.27e-02 1: 2.33e-02
4096
1.29e-03
1.26e-03 +: 2.95e-05
3.43e-03 3.47e-03
8192
5.66e-04
5.59e-04 2.86e-05
5.18--04
5.58e-044 6.19e-05
Low
32
2.38c-02
2.920-02
2.996-02:7..76-033
5,41c-02 1.16-02214446-00 1.46c-02
Volatility
64
6.57e-03
7.30e-00
2.7e-02244.733e
2.44e-02 8.32e-03
128
2.17e-03
3.34e-03
2.98c-03
+
1.10e-02 4.18e-03
7.56e-04
1.32e-03 1.51e-03
1.41e-03
--
4.81e-03 + 1.72e-03
512
2.85e-04
4.95e-04
5.44e-04
3.59c-04 1.07e-04
1.5le-03 + 6.71e-04
1024
1.13e-04
1.99e-04
= 0.66e-00
1.32e-04 1.22-04:..122 44.12e-05
5,60e-04 1 2.97e-04
2048
5.21e-05
4.87e-05
= 1.38e-06
5.03e-05 + 2.17e-06
1.92e-04 +: 1.85e-04
4096
2.32c-05
2.22c-05 7.64c-07
2.23e-05 8.73e-07
3.53e-05
+ 2.18c-05
8192
1.01c-05
9.66e-06 + 3.41c-07
9.40c-06 + 4.93c-07
1,00c-05 + 6.33c-07
Table 10: The table compares the standard Sobol' sequence (Joe & Kuo) against three variants: our implementation (LLM), our implementation with the flash-lite model (LLM-lite), and our im- plementation starting from random initialization (LLM random init.) across all six Asian option scenarios. We report the Mean Squared Error (MSE).
25

---

Published as a conference paper at ICLR 2026
Scenario
So
K
Suree
Training Example
50.00
45.00
0.3
7.06
Out-of-the-Monee
50.00
60.00
0.3
1.02
Al-the-Money
50.00
52.50
0.3
2.98
In-the-Money
50.00
40.00
0.3
11.02
High Volatility
50.00
52.50
6.43
Low Volatility
50.00
52.50
0.1
0.69
Table 11: Asian option scenarios used for testing. The Training Example used in the evaluation routine of the evolutionary search, All options have = 1.0 and r = 0.05.
Initial Price (So)
Strike Price (K)
Option Type
Scenario
Volatility (a)
Other Parameters
Asian
Training Example
50.00
45.00
0.3
Out-of-the-Monee
50.00
60.00
0.3
At-the-Money
50.00
52.50
0.3
In-the-Money
50.00
40.00
0.3
High Volatility
50.00
52.50
0.6
Volatility
50.00
52.50
0.1
Lookback
Base
100.00
0.2
High Volatility
100.00
0.4
Barrier
Base
100.00
100.00
0.2
Barrier Level: 85.00
Close Barrier
100.00
100.00
0.2
Barrier Level: 95.00
Basket (32D)
Low Correlation
100.00
100.00
0.2
DE 0.1
High Correlation
100.00
100.00
0.2
P: 0.8
Mixed Volatility
100.00
100.00
UO0.15,0.0 0.4)
P: 0.5
to-ofttte-Mooeyy
100.00
110.00
0.2
pr
Bermudan
At-thc- Miche-Money
100.00
100.00
0.2
Exercise Dates:
In-the-Monee
90.00
100.00
0.2
Exercise Dates:
Configuration Parameters for All Tested Option Scenarios. This table details the parame- ters for the 32-dimensional options used in the primary benchmark and generalizability tests. For all scenarios, the time to expiration (T) is 1.0 year and the risk-free interest rate (r) 0.05, An em-dash (-) indicates parameter is not applicable to that option type.
Table 12:
26

---

Published as a conference paper at ICLR 2026
Scenario
N
Squared Bias
Variance
MSE
p-value
Sobol
LLM
Sobol
LLM
Sobol
LLM
Training
32
2.09e-05
5,59c-07
0.2484
0.2523
0.2484
0.2523
0,9757
Example
64
5.37c-06
9.58e-06
0.0642
0.0665
0.0642
0.0665
0,9980
128
5.10e-06
4.96c-06
0.0183
0.0192
0.0183
0.0192
0.9999
256
2.87e-07
1.43e-07
0.0056
0.0061
0.0056
0.0061
1,0000
512
6.09e-08
3.51e-09
0.001646
0.001614
0.001646
0.001614
0.2841
1024
7.89e-08
1.33e-08
0.000542
0.000524
0.000542
0.000524
0,0548
2048
2.60e-08
2.30e-12
0.000233
0.000225
0.000233
0.000225
0.0199
4096
9.11e-10
3.46c-09
9.876c-05
9.346e-05
9.876c-05
9.346e-05
0.0058
8192
1.55c-09
1.93e-10
4.523c-05
4.104e-05
4.523c-05
4.104e005
7.0e-06
Out-of-the
32
3.73e-05
6.41e-07
0.1920
0.1979
0.1920
0.1979
0.9757
Money
64
1,48c-06
9.78e-09
0.0605
0.0631
0.0605
0.0631
0,9980
128
3.74c-07
1.66e-06
0.0220
0.0231
0.0220
0.0231
0.9999
256
3.07e-07
2.60e-07
0.0086
0.0089
0.0086
0.0089
1.0000
512
4.47c-08
4.05e-11
0.003469
0.003311
0.003469
0.003311
0.0208
1024
1.29e-08
4.45c-08
0.001457
0.001313
0.001457
0.001313
3,5e-08
2048
2.92c-11
3.36e-08
0.000619
0.000535
0.000619
0.000535
1,5c-14
4096
5.09e-08
1.17e-09
0.000258
0.000243
0.000258
0.000243
4.5c-04
8192
3.71e-09
3.24c-09
0.000117
0.000107
0.000117
0.000107
1.0e-06
At-the Money
32
5,90e-05
1.35e-09
0.3246
0.3346
0.3246
0.3346
0.9757
64
1.35e-05
1.17e-05
0.0873
0.0908
0.0873
0.0908
0,9980
128
9.72e-06
1.45-006
0.0277
0.0282
0.0277
0.0282
0.9999
256
1.12c-06
2.06e-08
0.0092
0.0098
0.0092
0.0098
1.0000
512
1,83e-08
5.57e-10
0.003248
0.003106
0.003248
0.003106
0.0208
1024
3.60c-08
4.58e-09
0.001228
0.001 168
0.001228
0.001168
0.0304
2048
7.85c-11
1.82e-08
0.000550
0.000521
0.000550
0.000521
0.0131
4096
3.16e-10
8.98e-10
0.000240
0.000226
0,000240
0.000226
0.0019
8192
7.47e-09
3.31c-09
0.000102
9.522c-05
0.000102
9.523c-05
0.0100
In-the
32
1.93e-05
3.11e-09
0.1397
0.1428
0.1398
0.1428
0.9757
Money
64
2.30e-06
2.86c-06
0.0367
0.0382
0.0367
0.0382
0,9980
128
1.67e-06
1.74c-06
0.0106
0.0111
0.0106
0.0111
0,9999
256
1.08c-09
3.54c-08
0.003083
0.003373
0,003083
0.003373
1,0000
512
3.07e-10
4.26c-08
0.000868
0.000859
0.000868
0.000859
0.2841
1024
1.53e-09
6.44e-10
0.000261
0.000260
0.000261
0.000260
0,2581
2048
3.54e-10
1.28c-08
0.000101
0.000100
0.000101
0.000100
0.3878
4096
2.15e-11
3.38e-09
4.037e-05
3.949c-05
4.037e-05
3.949c-05
0.2180
8192
2.77e-10
3.70e-09
1.766e-05
1.682c-05
1.766e-05
1.683e005
0.0047
High
32
0.000333
2.83e-07
2,0813
2.1493
2,0816
2.1493
0.9757
Volatility
64
5,87c-05
4.32e-05
0.5713
0.5958
0.5714
0.5958
0,9980
128
4.41c-05
5.25e-06
0.1810
0.1877
0.1811
0.1877
0.9999
256
1.78e-06
1.90c-06
0.0598
0.0643
0.0598
0.0643
1.0000
512
2.34c-09
1.49e-08
0.0205
0.0196
0.0205
0.0073
1024
1.21c-07
2.19e-07
0.007565
0.007060
0.007565
0.007060
3,4c-04
2048
6.34c-08
1.22e-07
0.003145
0.002867
0.003145
0.002867
1.8c-06
4096
2.93e-10
6.09e-09
0,001293
0.001233
0,001293
0.001233
0.0022
8192
3.59e-08
4.35e-08
0,000566
0.000532
0.000566
0.000532
0.0156
Low
32
5.75e-06
1.26e-08
0.0238
0.0246
0.0238
0.0246
0.9757
Volatility
64
1.28c-06
6.48e-07
0.0066
0.0068
0.0066
0.0068
0,9980
128
7.71c-07
2.71e-08
0.002173
0.002205
0.002174
0.002205
0.9999
256
1.08e-07
1.12e-08
0.000756
0.000803
0.000756
0.000803
1,0000
512
6.51e-09
1.55e-09
0.000285
0.000275
0.000285
0.000275
0.0208
1024
5.18e-09
1.20e-09
0.000113
0.000109
0.000113
0.000109
0.0717
2048
2.66e-10
2.17e-09
5.214e-05
4.961c-05
5.214e-05
4.961c-05
0,0056
4096
1.26e-10
1.39c-10
2.323c-05
2.196e-05
2.323c-05
2.196e-05
5.7e-04
8192
3.43e-10
3.09e-10
1. 1.008e-05
9.300e-06
1.008e-05
9.300e-06
1.3e-04
Table 13: Integration results for the Asian option across all six scenarios. The table compares the standard Sobol' sequence (Joe & Kuo) against the direction numbers discovered by L.M evolu- tionary search. We report the Mean Squared Error (MSE), its constituent parts (Squared Bias and Variance), and the FDR-corrected p-value from a one-sided Wilcoxon signed-rank test. P-values below 0.05 and the best methods bolded.
27

---

Published as a conference paper at ICLR 2026
Scenario
Squared Bias
Variance
MSE
p-value
Sobol
L.M
Sobol
L.M
Sobol
LLM
Base
32
3.90e-06
8.56e-05
1.2138
1.2443
1.2138
1.2444
0.9968
64
4.16e-06
3.01e-05
0.4725
0.4839
0.4725
0.4840
1.0000
128
3.146-06
2.35e-07
0.1272
0.1352
0.1272
0.1352
1.0000
256
1.6le-06
4.90c-06
0.0352
0.0381
0.0352
0.0381
1.0000
512
8.86e-07
1.10c-06
0.0121
0.0120
0.0121
0.0120
0,4675
1024
9.64c-07
1.43e-06
0.005021
0.004905
0.005022
0.004907
0.0342
2048
6.88c-08
1.29c-07
0.001834
0.001744
0.001834
0.001744
4.1c-04
4096
1.32e-07
5.59c-08
0.000807
0.000768
0.000807
0.000768
2.4c-04
8192
5.21e-08
3.61e-08
0.000407
0.000377
0.000407
0.000377
9.8e-09
High
32
2.13e-04
4.44c-04
8,5591
8.7869
8,5593
8.7873
0.9968
Volatility
64
6.04e-08
1.16e-04
3.3964
3.5134
3.3964
3.5136
1.0000
128
2.87e-05
2.57e-06
0.8854
0,9982
0.8854
0,9982
1,0000
256
3.14e-06
1.46c-05
0.2403
0.2666
0.2403
0.2663
1,0000
512
7.96c-07
1.11e-07
0.0809
0.0792
0.0809
0.0792
0,3836
1024
2.57e-06
7.04c-06
0.0325
0.0312
0.0325
0.0312
6.2e-04
2048
3.37e-07
4.92e-08
0,012003
0.011234
0.012003
0.011234
2.3e-06
4096
1.47e-07
2.32e-08
0,004862
0.004558
0.004862
0.004558
2.7e-06
8192
2.12e-10
6.25e-09
0.002335
0.002057
0.002335
0.002057
9.8c-17
Table 14: Integration results for the Lookback Option across two scenarios. We report FDR-
corrected P-values from a one-sided Wilcoxon signed-rank test. P-values below 0.05 and the best methods are bolded.
Scenario
Squared Bias
Variance
MSE
p-value
Sobol
LLM
Sobol
LIM
Sobol
LLM
Base
32
1.13e-04
2.87e-05
2.2540
2.3339
2.2541
2.3339
0.9968
64
2.17e-05
4.56e-06
0.8386
0.8894
0.8386
0.8894
1.0000
128
2.32e-05
2.37e-05
0.2533
0.3101
0.2534
0.3101
1,0000
1,97e-06
7.40e-06
0.0707
0.0834
0.0707
0.0834
1,0000
512
5.72e-07
4.84c-06
0.0245
0.0241
0.0245
0.0241
0.4675
1024
7.63e-07
7.40e-07
0.009948
0.009819
0.009948
0.009819
0.3313
2048
4.03e-07
1.12e-06
0.004277
0.004230
0.004277
0.004231
0.3025
4096
8.61e-08
1.00e-07
0.002019
0.001992
0.002019
0.001992
0.3863
8192
1.82e-08
5.70e-08
0.000985
0,000955
0.000985
0.000955
0.0859
Close
32
9.55c-04
1.27e-04
3,2150
3.2754
3.2160
3.2755
0,9968
Barrier
64
6.55c-05
1.44e-05
1.1229
1.2054
1.1229
1.2054
1.0000
128
5.52c-05
1.14e-05
0.4541
0,5183
0.4542
0.5183
1,0000
256
3.62e-05
8.39c-07
0.1733
0.1775
0.1733
0.1775
1,0000
512
3.76e-05
7.16c-07
0.0723
0.0690
0.0723
0.0690
0.2364
1024
2.27e-05
2.62e-06
0.0328
0.0326
0.0328
0.0326
0.4459
2048
1.52e-06
1.13e-09
0.01594
0.01648
0.01594
0.01648
0.9126
4096
3.83e-07
4.25e-07
0.007027
0.006727
0.007028
0.006728
0.0574
8192
3.48e-09
3.79e-07
0.003379
0.003266
0.003379
0.003267
0.1613
Table 15: Integration results for the Barrier Option across two scenarios. We report FDR-corrected p-values from a one-sided Wilcoxon signed-rank test. P-values below 0.05 and the best methods are bolded.
28

---

Published as a conference paper at ICLR 2026
Scenario
Squared
Bias
Variance
MSE
p-value
Sobol
LLM
Sobol
LLM
Sobol
LLM
Low
32
3.82e-06
5.13e-07
0.115
0.118
0.115
0.118
0.997
Correlation
64
6.60e-07
1.99e-06
0.0342
0.0358
0.0342
0.0358
1.00
128
1.18e-06
1.38e-06
0.0109
0.0126
0.0109
0.0126
1.00
256
5.16e-08
1.38e-08
0.00360
0.00390
0.00360
0.00390
1.00
512
3.04c-07
1.15e-08
0.00126
0.00123
0.00126
0.00123
0.238
1024
2.00c-07
1.79e-08
5.04c-04
4.73e-04
5.04c-04
4.73e-04
0.0051
2048
5.00c-08
6.46e-11
2.29c-04
2.14e-04
2.29c-04
2.14e-04
3.8e-04
4096
1.12e-10
1.32e-09
9.89e-05
9.04c-05
9,89e-05
9.04c-05
1.7c-04
8192
1.39e-09
1.72c-10
4.74e-05
4.16e-05
4.74e-05
4.16e-05
1.9c-11
High
32
9.26c-05
9.80e-05
0.248
0,256
0,248
0.256
0.997
Correlation
64
7.77e-06
1.10c-05
0.0602
0.0570
0.0602
0.0571
2.6e-05
128
6.94c-06
2.65e-06
0.0163
0.0168
0.0163
0.0168
1.00
256
7.96e-07
8,54c-07
0.00491
0.00494
0.00491
0.00495
1.00
512
6.53e-08
8.95e-08
0.00153
0.00152
0.00153
0.00152
0.468
1024
3.80e-10
1.16e-08
4.63e-04
4.47e-04
4.63e-04
4.47e-04
0.0784
2048
3.45e-09
1,57e-09
1.75e-04
1.78e-04
1.75e-04
1.78c-04
0.742
4096
8.68e-10
7.39e-09
7.13e-05
6.61e-05
7.13e-05
6.61e-05
8.2e-04
8192
6.90e-09
1.17e-10
2.60e-05
2.45e-05
2.60e-05
2.45e-05
0.00530
Mixed
32
2.19e-04
9.65e-05
0.775
0.799
0.776
0.799
0.997
Volatility
64
9.78e-06
1.28c-05
0.210
0.207
0.210
0.207
0.540
128
1.95c-05
2.19e-06
0.0647
0.0702
0.0647
0.0702
1.00
256
4.42c-06
1.52e-06
0.0211
0.0216
0.0211
0.0216
1.00
512
6.51e-07
6.07c-07
0.00620
0.00611
0.00620
0.00611
0.640
1024
1.17e-08
5,20c-09
0.00209
0.00203
0.00209
0.00203
0.0784
2048
1.09e-08
4.65e-11
9.47e-04
9.32e-04
9.47e-04
9.32c-04
0.0672
4096
3,40e-09
2.79c-08
4.47e-04
3.93e-04
4.47e-04
3.93e-04
1.5e-10
8192
1,80c-08
2.19e-12
1.56e-04
1.43e-04
1.56c-04
1.43e044
1,5e-05
Out-of-the
32
2.43e-05
2.29e-09
0.117
0.121
0.117
0.121
0.997
Money
64
2.81e-06
4.90c-09
0.0367
0,0387
0.0367
0.0387
1.00
128
2.55e-06
2.01e-07
0.0127
0.0143
0.0127
0.0143
1,00
256
5.97c-07
4.26e-10
0.00437
0.00463
0.00437
0.00463
1.00
512
2.65e-07
2.67e-08
0.00160
0.00155
0.00160
0.00155
0.236
1024
3.70e-08
2.12e-09
6.52e-04
6.23e-04
6.52e-04
6.23e-04
0.0651
2048
5.68e-09
3.00e-11
2.92e-04
2.81e-04
2.92e-04
2.81c-04
0.0856
4096
1.01e-08
5.67e-09
1.25e-04
1.20e-04
1.25e-04
1.20e-04
0.0027
8192
1.58e-11
1.66c-09
5.95c-05
5.39e-05
5.95c-05
5.39e-05
2.5e-06
Table 16: Integration results for the 32-dimensional Basket Option across four scenarios. We report FDR-corrected p-values from a one-sided Wilcoxon signed-rank test. P-values below 0.05 and the best methods are bolded.
29

---

Published as a conference paper at ICLR 2026
Scenario
N
Squared Bias
Variance
MSE
p-value
Sobol
LLM
Sobol
LLM
Sobol
LLM
At-the
32
0.817
0.827
0.574
0.583
1.39
1.41
0.997
Money
64
0.254
0.263
0.240
0.237
0.495
0,500
1.00
128
0.0698
0.0725
0.0941
0.0951
0.164
0.168
1.00
256
0.0182
0.0192
0.0441
0.0445
0.0624
0.0637
1.00
512
0,00440
0.00410
0.0218
0.0221
0.0262
0.0262
0,468
1024
9.87e-04
9.33e-04
0.0118
0.0112
0.0128
0.0122
0.0115
2048
2.77e-04
2.46e-04
0.00596
0.00571
0.00623
0.00596
0.0018
4096
5,23e-05
4.25e-05
0.00270
0.00268
0.00275
0.00272
0.386
8192
1.10e-05
8.45e-06
0,00132
0.00127
0.00134
0,00128
0.0457
In-the
32
0.889
0.898
0.563
0.555
1.45
1.45
0.997
Money
64
0.223
0.230
0.251
0.251
0.473
0.481
1.00
128
0.0478
0.0489
0.118
0.117
0.166
0.166
1.00
256
0.00920
0.00930
0.0564
0.0563
0.0656
0.0656
1.00
512
0.00200
0.00190
0.0284
0.0278
0.0304
0.0297
0.468
1024
3.75e-04
4.09c-04
0.0145
0.0135
0.0149
0.0139
0.0142
2048
6.75e-05
7.12c-05
0.00735
0.00658
0.00742
0.00666
1.2e-05
4096
1.17e-05
9.36e-06
0.00337
0.00313
0,00338
0.00314
7.4e-04
8192
3.56e-06
3.86e-06
0.00166
0.00152
0.00166
0.00153
1.6e-04
Table 17: Integration results for the Bermudan Option corrected p-values from a one-sided Wilcoxon signed-rank test. P-values below 0.05 and the best methods are bolded.
scenarios. We report FDR-
30


Reviews:
[{"review_id": "KoRVmnFex7", "reviewer": "Reviewer_2Cbf", "summary": "The paper applies OpenEvolve, a Large Language Model (LLM)-guided evolutionary program synthesis based on AlphaEvolve, to improve low-discrepancy sequences. Specifically, they apply this methodology to (1) construct finite 2D and 3D point sets with minimal star discrepancy (leading to small integration error); (2) optimising the direction numbers of Sobol’s sequence to reduce high-dimensional integration error \nThe LLM-guided evolutionary approach generates and mutates code that produces candidate point sets or sequence parameters, guided by a fitness function that measures, respectively, star discrepancy or integration error. \nSpecifically, in 2D, for a fixed number of points $N \\leq 10$ the method recovered previously known optimal point sets and found new point sets with lower star discrepancy than previously discovered for $N > 30$ (0.0150 at $N=100$ vs the prior 0.0188).\nFor 3D sequences, the method matched all known optima up to N=7, and found new record-low discrepancies for N>8, for which optima are unknown. \nFinally, for Sobol sequences in 32 dimensions, the evolved direction parameters reduced the mean squared error of randomised QMC integration on an Asian option pricing task compared to the standard Joe-Kuo parameters.", "questions": "* Justify the choice of OpenEvolve instead of other LLM-guided evolutionary algorithms\n\n*  Which LLM was used, what inference settings, and why?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "* The paper applies a powerful emerging method (LLM-guided evolutionary program synthesis) to a well-suited, long-standing problem in a creative manner.  This is very significant in two ways. Firstly, it is an interesting application of LLM-guided evolutionary methods, which are establishing themselves as an important tool in scientific discovery. Secondly, discovering new low-discrepancy point sets and sequence parameters is significant for the QMC and numerical methods community.  \n\n* The paper is clearly organised, provides an adequate amount of background, and contains extensive experimental results to back up its claims. The authors compare their results against a thorough set of baselines, including traditional sequences (Halton, Sobol, etc), simple lattice heuristics, and even more recent GNN-based MPMC methods. \nWhile the Sobol sequence parameters were optimised with respect to the integration error of Asian option pricing, the authors evaluate their performance on multiple other exotic option types, showing that the solution generalises beyond this specific problem.", "weaknesses": "* Presentation:  line 267 has a missing reference. In addition, line 097 would benefit from a citation when referring to the “total variation of $f$ in the sense of Hardy and Krause”. \n\n* The paper would benefit from more details on the LLM settings (which LLM is used?) and, perhaps, comparisons across multiple LLMs or LLM ensembles. I would expect that different LLMs would lead to different programs and solutions of different quality/optimality.\n\n* A comparison or discussion justifying the choice of OpenEvolve over other LLM-guided evolutionary methods (such as ShinkaEvolve) would be interesting. \n\n* Finally, one limitation compared to the previously established baseline Sobol parameters by Joe-Kuo is the solution’s potential specialisation for the options pricing problem.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "PiAbUOey8B", "reviewer": "Reviewer_ihGh", "summary": "This paper applies an LLM-guided evolutionary program synthesis framework (OpenEvolve) to two problems in Quasi-Monte Carlo (QMC) design: \n\n**1**. constructing finite 2D/3D point sets with low star discrepancy \n\n\n**2**. optimizing Sobol' direction numbers to reduce integration error for high-dimensional financial models. \n\nThe authors report finding new state-of-the-art 2D point sets for $N \\ge 40$ and discovering new Sobol' parameters that outperform the standard Joe-Kuo parameters on a suite of financial option pricing tasks.", "questions": "I have two main questions, which are based on the weaknesses above:\n\n**1.** Corresponding to Major Weakness 1, how could the authors justify that the improvement in the point-set problem comes from the LLM-guided approach, rather than just the application of scipy.optimize?\n\n**2.** For the Sobol' discovery, how much credit belongs to the LLM's \"intelligent search\" versus the fact that it started with the Joe & Kuo parameters as its initialization? Have the authors tried running the search from a random initialization to see if it can discover good parameters from scratch?", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 1, "strengths": "**1**. The paper successfully refined the Joe&Kuo's result, finding new Sobol's direction numbers that result in a statistically significant reduction in integration error (MSE) for a suite of 32-dimensional financial integration tasks;\n\n**2**. The paper introduces the current LLM tools for solving a long-standing discrete optimization problem;\n\n**3**. The paper is well-written, clearly organized and does a good job of introducing the complex technical details of QMC.", "weaknesses": "**Major Weaknesses:**\n\n\n**1**.Minimal Contribution in Point Set Discovery:  The paper's \"two-phase\" strategy for point sets is a significant weakness. As shown in Figure 1, the LLM's \"direct construction\" in Phase 1 provides almost no improvement (0.0962 $\\rightarrow$ 0.0924). The entire significant gain (0.0924 $\\rightarrow$ 0.0744) comes from Phase 2, which is just the LLM generating code to call a standard classical optimizer (scipy.optimize.minimize). If this is the case, then it is not a novel discovery by the LLM; it's the automation of a standard workflow any human researcher would perform.\n\n**2**. Lack of Methodological Novelty: The core evolutionary framework is a direct application of the pre-existing OpenEvolve. The paper does not propose any new methodological innovations to this framework. Its contribution is limited to applying this existing tool to the QMC domain.\n\n**3**. Highly Incremental Sobol' Results: While statistically significant, the improvement in the Sobol' task is extremely small in absolute terms (e.g., an MSE of 4.52e-05 vs. 4.10e-05). Crucially, the search was initialized with the strong Joe & Kuo (2008) parameters. This frames the discovery not as a major breakthrough, but as a minor, local refinement of an existing solution.\n\n\n**Minor Weakness:**\n\n**1**. The authors acknowledge that they performed only one evolutionary run per problem. Since evolutionary algorithms are stochastic, this single run provides good new parameters but tells us nothing about the robustness or reliability of the search method.\n\n\n\nIf there's any misunderstanding, I would be more than happy to correct my opinions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Axb12BcSGx", "reviewer": "Reviewer_qCAy", "summary": "The paper applies LLM to perform numerical integration with an evolutionary computation angle to improve the integration results.\nOverall a nice paper.", "questions": "- P values to the order of 10^-14 indicate a  population with very little variance, is there any reason for such lack of variance?\n\n- I think, this is a very interesting application, However, one needs to understand and study the computational complexity of getting such solutions. Moreover, I get that there are about 3 applications studied in this paper, I would like to know that generality of this approach, do you have to fine tune the LLM in some way.\n\n- However, about the mutation, does the LLM involve some kind of RL based training to generate better mutations.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "Well written paper, \nThe application of LLM with an evolutionary loop to generate populations is interesting\nThe method has an evolving loop of introducing mutations to provide better solutions.", "weaknesses": "It appears that the approach provides better solution with N becomes larger and larger. However, the improvement is marginal \nOver SOBOL and clement et. Al in the sense of MSE. While this marginal improvement is great, the computational complexity of \nRunning an evolutionary loop is something unknown right now.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "MMfyNeo4fZ", "reviewer": "Reviewer_W9ms", "summary": "This paper introduces an LLM-guided evolutionary framework to automate the discovery of quasi-Monte Carlo (QMC) constructions. The authors treat the design of low-discrepancy point sets and Sobol’ sequence direction numbers as a program synthesis problem. Within the proposed OpenEvolve framework, an LLM iteratively mutates Python programs that generate candidate point sets or Sobol’ parameters, guided by fitness scores (e.g., inverse star discrepancy or rQMC mean squared error).", "questions": "See Weakness.", "limitations": "", "rating": 4, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. Casting QMC design as program synthesis is conceptually elegant and connects symbolic LLM reasoning with continuous numerical optimization.\n\n2. The paper carefully combines constructive heuristics with iterative optimization, and the experimental evaluation uses statistically robust paired tests", "weaknesses": "1. The “LLM-guided” aspect is somewhat opaque. How much improvement stems from the LLM’s structured code editing versus brute-force evolutionary search or the built-in optimizer (SLSQP)?\n\n2. This problem setting seems too narrow, and i am not an expert in this domain so i do not understand the significance of this improvement. \n\n3. Are these compared baselines the best baseline in this field? Is there some other learning-based baselines? '\n\n\n4. In my view, it seems the authors only carefully design a prompting procedure to solve a specfic problem. Not sure these design strategies can be used for other problems. \n\n5. The scale N/d is kind of limited, and we do not know whether this method can scale. In my understanding, LLM can be quite hard to scale for high-dimensional numbers.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "OMRNvg1Mk0", "reviewer": "Reviewer_jdvu", "summary": "This paper proposes a new methodology to find low-discrepancy point sets, by using an evolutionary algorithm where the mutation phase is done by an LLM. Two settings are inspected: in the first setting, the goal is to minimize the star discrepancy of a point set in 2D and 3D. Here the generation space is the space of Python programs. In the second setting, the goal is to minimize the integration error for a 32D option-pricing task. Here the generation space is the set of parameters (so-called direction numbers) for Sobol digital nets. In both cases, the solution found by the proposed method outperforms baselines.", "questions": "How many different direction numbers were evaluated when running the evolutionary algorithm on the option pricing problem? Is it 2000 as reported in Appendix B.3?", "limitations": "", "rating": 2, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "Guiding program synthesis through LLMs is a thriving avenue of research, which has shown progress in various fields. I appreciate that this paper brings this idea to the problem of low-discrepancy point set generation. The experiments show an improvement over the best known star-discrepancy sets in 2D and 3D, and improve over Joe and Kuo direction numbers in the Asian option pricing experiment.", "weaknesses": "My main concern is the lack of scientific evaluation, especially given that the proposed method is highly resource-intensive, involving thousands of LLM calls, each requiring computing the fitness of the point set. Such evaluation is critical to assess the proposed approach. Examples of lacking scientific evaluation include:\n- comparing with simpler methods for exploring the space of low-discrepancy sequences. I discuss this point below for each experimental setting.\n- reporting the improvement in terms of fitness during the run of the evolutionary algorithm.\n- comparing with simple (single- or multi-turn) LLM prompting without population-based evolution.\n\nFor the star-discrepancy experiment, a natural baseline is local optimization from an existing high-quality construction. Since the two-phase prompting strategy already instructs the LLM to “use scipy optimization routines such as scipy.optimize.minimize,” it would be appropriate to compare against standard optimizers (e.g., L-BFGS-B, SLSQP) initialized from Clément et al. or other RQMC sequences, using a comparable computational budget (e.g., 2000 randomizations).\n\nFor the Asian call experiment, the LLM modifies direction numbers in dimensions 4–6 relative to the Joe and Kuo numbers. Since the search space over these dimensions has size 2^12, comparable to the 2000 LLM calls reported, it is unclear if the proposed algorithm provides an advantage over random search. Including such comparison (and comparing to exhaustive search to get an indication of how far away the solution found by the LLM is from the optimum) would strengthen the analysis. Another baseline is the LatNetBuilder software, which is specifically designed for identifying good low-discrepancy point sets, including Sobol sequences. An interesting approach would be to rank direction numbers in terms of some fitness measure (e.g., star-discrepancy t-value) using the software, then to evaluate them in this order on the Asian call experiment, and to compare the performance with the LLM-guided evolutionary approach, as a function of the number of tested direction numbers.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)