Decision: the initial design of Bayesian optimization

[tachsin]

Batch B of docs/optimization-plan.md leaves the default initial design size of Bayesian optimization (0.13) to settle.

Option	Source	For n = 6
10n	Jones, Schonlau and Welch (1998); Loeppky, Sacks and Welch (2009)	60
2(n + 1)	Common in BO libraries	14

With a budget of 100 evaluations in 6 dimensions, 10n spends 60% of it before the model guides anything. With 2(n + 1), the first hyperparameter fits use only 14 points.

Recommendation: 2(n + 1), with a setting to change it, and the docs citing Loeppky et al.'s 10n for when the budget allows. The examples report how they reach the optimum with the default.

Also to decide: whether the design is a Latin hypercube only (McKay et al. 1979, added in the 0.13 draft PR), or offers a scrambled Sobol sequence too (Joe and Kuo 2008). Recommendation: Latin hypercube now, Sobol when someone needs it.

[tachsin]

I tested this with a prototype written from the textbook definitions (GPML and Ament et al.), in numpy and scipy:

• a GP with a constant mean and an ARD Matérn 5/2 kernel, hyperparameters by maximum likelihood with L-BFGS-B, refitted every iteration;
• log-EI maximized from many starts;
• a Latin hypercube design.

Its gradients are checked against finite differences. I ran 30 seeds per configuration and 100 evaluations per run, 1,320 runs in all.

Simple regret (median) and evaluations to reach f + 1e-3* (p: Mann-Whitney against 2(n+1) at 100 evaluations, uncorrected):

problem	n+1	2(n+1)	5n	10n
Branin, evaluations to target (all 30 reach it)	28	29	30	35
Six-hump camel, regret at 100	7.2e-5 (p 0.004)	1.9e-5	2.4e-5	1.8e-5
Six-hump camel, evaluations to target	49	49	52	57
Goldstein-Price (log-transformed for the model), evaluations to target	36	33	31	41
Hartmann 3, runs reaching the target / median evaluations	19 / 35	21 / 32	22 / 27	27 / 36
Hartmann 6, regret at 50 / at 80	6.2e-2 / 8.0e-5	1.1e-3 / 6.7e-5	1.2e-1 / 7.1e-5	(design only) / 4.8e-3
Hartmann 6, evaluations to target	56	48	58	82
Ackley 5-D, regret at 50 / at 100	11.1 / 2.49 (p 0.010)	5.5 / 2.19	7.6 / 2.19	17.8 / 2.05
Rosenbrock 5-D, regret at 100	302	180	337	1260 (p < 0.001)

• 2(n+1) is best or tied everywhere, and never significantly worse than another size at 100 evaluations.
• 10n wastes budget: 7× the regret on Rosenbrock 5-D, and the most evaluations to the target (82 against 48 on Hartmann 6).
• n+1 is the risky minimum: significantly worse on Six-hump and Ackley.
• 5n is indistinguishable at 100 evaluations, but worse early.

The first model is poor with small designs, but that doesn't hurt the search. With n+1 or 2(n+1) points in 3 or more dimensions, most first fits push length scales to their upper bound: a near-constant model, whose predictions are no better than a constant (normalized RMSE ≈ 1). Refitting each iteration recovers it, as the table shows.

The output transform matters more than the design size. On raw Goldstein-Price, no design size reached the target within 100 evaluations; with log-transformed outputs, all did in about 33. The plan has no output warping. A log or rank option, or at least a docs note, belongs in batch B.

Latin hypercube against scrambled Sobol against random: no significant difference (all p ≥ 0.17), agreeing with Bossek et al. below.

Literature (texts read unless marked):

• Jones, Schonlau and Welch (1998), §4: "about n = 10k points are needed in the initial design".
• Loeppky, Sacks and Welch (2009): the 10d rule is about the GP's prediction accuracy, not optimization. In their words, "optimization of a target criterion could raise other considerations".
• Bossek, Doerr and Kerschke (GECCO 2020): "a clear tendency towards small initial ratios"; design types matter little.
• Brockhoff, Bischl and Wagner (2015, abstract): "large initial samples seem detrimental".

Correction to the issue: "2(n+1) is common in BO libraries" overstates it. Defaults range from 4 to 10 points, d+1 (HEBO), 2d (Ax, TuRBO's recommendation), 4d (mlrMBO) and 5d (Dragonfly). 2(d+1) is MATSuMoTo's default and is used in pySOT's and BoTorch's examples.

Recommendation (confidence medium-high for up to 6 dimensions and budgets of 40-100; I didn't test more than 10 dimensions):

• 2(n+1) by default, as a setting.
• The docs cite Loeppky et al.'s 10d for when an accurate global model is the goal.
• Latin hypercube now, Sobol only on request.
• An output transform option goes into batch B.

[tachsin]

Settled for batch B (0.13), following the evidence above:

• Initial design: 2(n + 1) points by default, as a setting (initial_points). The docs cite Loeppky et al.'s 10n for when an accurate global model is the goal.
• Design: a Latin hypercube (McKay et al. 1979). A scrambled Sobol sequence comes only on request.
• Output transform: an option in batch B, for objectives such as Goldstein-Price that the model fits better on a log scale.
• Names (deferred here from Decision: names and defaults of the local and Bayesian methods #376): Bo, like Ga, De and Es; the public Gaussian process in model::gp; and the acquisition functions in algorithm::bo::acquisition, with Bo, as the draft PR feat(bo): Bayesian optimization with a Gaussian process model, log-EI and a Latin hypercube design #382 has them.