Example 11: expected improvement with SMC search

Script: examples/example11_optim_EI_smc.py

Purpose

The script runs the same one-dimensional EI loop as example10, but replaces the fixed candidate grid by an SMC particle set. The particle set is moved toward regions where the posterior probability of improvement is not negligible. This follows the general idea of using SMC particles to represent a sequence of search targets, as in Bayesian subset simulation [1], while using the EI criterion introduced for efficient global optimization [4].

What is computed

• posterior mean and variance at the current particle positions.

• EI values on the particle set.

• a target log-density proportional to log P(-Y(x) > -min(zi) | observations) inside the input box.

• SMC particle updates using reweighting, resampling, and Markov moves.

• one new objective evaluation per sequential step.

Main objects

• gpmpcontrib.optim.expectedimprovement.ExpectedImprovementSMC

• gpmpcontrib.SequentialStrategySMC

• gpmp.mcmc.smc.SMC

Outputs

Run python examples/example11_optim_EI_smc.py from the repository root to execute the example. Regenerate the static figure with cd docs && python make_example_results.py.

Top panel: GP posterior from the current observations. Middle panel: EI on a dense plotting grid. Lower panel: improvement-probability profile with SMC particles. Particle heights show the target density up to a multiplicative constant. Particles concentrate where the improvement probability is not close to zero.

Source excerpt

# -- define a mono-objective problem

problem = gpc.ComputerExperiment(
    1,  # dim of search space
    [[-1], [1]],  # search box
    single_function=gp.misc.testfunctions.twobumps,  # test function
)


# -- create initial dataset

nt = 2000
xt = gp.misc.designs.regulargrid(problem.input_dim, nt, problem.input_box)
zt = problem(xt)

ni = 3
ind = [100, 1000, 1600]
xi = xt[ind]

# -- initialize a model and the ei algorithm
model = gpc.Model_ConstantMean_Maternp_REML(
    "GP1d",
    output_dim=problem.output_dim,
    mean_specification={"type": "constant"},
    covariance_specification={"p": 2},
)

eialgo = ei.ExpectedImprovementSMC(problem, model)

eialgo.set_initial_design(xi)