Scatter search gives better fits than differential evolution on our problems, and it is the method best suited to a cluster, because the work it does in one iteration is population_size * (population_size - 1) simulations that can all run at the same time. A reference set of 15 fills 210 processors without anyone having to invent anything.
This epic has one theme. Finish the scatter search template as Fred Glover described it, then make the method understand that a simulation of a stochastic model gives a noisy answer.
The four steps below can ship one at a time and in this order. Steps 1 and 2 improve the method for everyone, including people fitting deterministic models. Steps 3 and 4 are for stochastic models and are the part no other tool does.
Background on what we have today
The implementation is in pybnf/algorithms/optimizers/scatter_search.py. It keeps a small set of good parameter sets called the reference set. Each iteration it combines every ordered pair of reference members to produce candidate parameter sets, simulates all of them, and replaces a reference member when one of its children scores better. A member that fails to improve for several iterations is retired to an archive and replaced with a fresh random point.
Against Glover's original description we have the diverse starting population, the systematic combination of every pair, and the combination formula itself, which comes from Egea and colleagues. Two pieces are missing.
Step 1: add the improvement method
Glover's template applies a local search to the candidates produced by combination. Egea's version does the same. Ours does no local search at all and relies purely on recombination. This is the largest gap and the most likely source of better fits.
We already have the part needed. SimplexRunner in pybnf/algorithms/optimizers/simplex.py is a self-contained, picklable state machine that runs one Nelder-Mead simplex search, works in parameter set space, and knows nothing about the trajectory, backup, or the job scheduler. PowellRunner in pybnf/algorithms/optimizers/powell.py is the same shape. Both were built for the concurrent multi-start work and can be driven from inside scatter search.
Local search is expensive, so it should not run on every candidate. Egea's version uses filters to decide when a candidate is promising enough to be worth refining. We need something similar. On a cluster the refinements can run at the same time as each other, so the cost is wall clock only when we are short of processors.
Step 2: choose the diverse half of the reference set by distance
round_1_init at line 146 builds the first reference set by taking the best half by score and then filling the rest with rng.choice over the remainder. Glover's point is that the second half should be the most diverse members, chosen by their distance from what is already in the set. Picking at random gives diversity only on average, and in a problem with many parameters that is much weaker than choosing for it deliberately.
This is a small change to one method and should be measurable on its own.
Step 3: make the reference set aware of noise
When a model is stochastic, running it twice with the same parameters gives two different objective values. Scatter search currently stores one simulation's result as a member's score and treats it as fact.
The useful observation is that scatter search never needs the objective values to be precise. It only needs the ordering to be right. The step size in the combination formula comes from the gap in rank between the two members being combined:
beta = (abs(hi - pi) - 1) / (self.popsize - 2)
Sorting the reference set, deciding whether a child beat its parent, and counting whether a member is stuck are all ranking decisions too.
So the reference set should hold an estimate and a measure of uncertainty for each member rather than a single number, and extra simulations should be spent only where the ordering is genuinely in doubt. Two members whose uncertainty ranges overlap, and whose rank gap is driving a combination, are worth more simulations. Members that are clearly separated are worth none. How many repeats a parameter set receives then follows from how contested its position is, rather than from a setting the user has to choose.
This also fixes a real fault. The stuck counter at line 197 does not know about noise. A member that got a lucky simulation cannot be beaten by honest children, so its counter climbs until it is retired into the local_mins archive as a local minimum that it never was. The archive of best results quietly fills with artifacts of the noise.
Related work is in #659, which fixes the same kind of error in the final answer at the end of a run. That issue is smaller, it stands alone, and it will tell us how large this problem actually is on real fits before we commit to this step.
Step 4: stop waiting for the whole iteration, and give idle processors something useful to do
ScatterSearch sets waits_for_full_generation = True, so every iteration waits for its slowest simulation. This is worst for exactly the models we care about here, because the running time of a stochastic simulation depends on the trajectory it happens to take, so the spread in running times is wide.
Once step 3 exists, the fix for the idle processors is already in hand. When a processor becomes free and no new candidate is ready, give it another repeat of a parameter set whose rank is still uncertain. The idle processor problem and the noise problem then share one answer.
How this differs from saCeSS
Penas and colleagues published a parallel, cooperative, self-adapting scatter search in 2017, and our source already notes that we implemented the plain version and not theirs. It is worth being clear about why this epic is not a copy of it.
Their method is built for large deterministic models, where the objective is a fixed number. A processor can either evaluate a new candidate or run another independent copy of the search. With a stochastic model there is a third choice, which is to reduce the uncertainty about a candidate already seen, and their design has no way to express it because it has no notion of noise.
They also run many independent copies of the search, each doing its combinations one after another. We run a single search and evaluate all of its combinations at once. Both are reasonable uses of a large machine, and ours is already built.
Also worth doing
The size of the reference set decides how many simulations run at once, so the right value is close to the square root of the number of processors available. We already read the number of connected workers in Algorithm._report_parallelism. We could set a sensible reference set size automatically instead of asking users to work out population_size * (population_size - 1) against their allocation by hand, which is what docs/cluster.rst tells them to do today.
Scatter search gives better fits than differential evolution on our problems, and it is the method best suited to a cluster, because the work it does in one iteration is
population_size * (population_size - 1)simulations that can all run at the same time. A reference set of 15 fills 210 processors without anyone having to invent anything.This epic has one theme. Finish the scatter search template as Fred Glover described it, then make the method understand that a simulation of a stochastic model gives a noisy answer.
The four steps below can ship one at a time and in this order. Steps 1 and 2 improve the method for everyone, including people fitting deterministic models. Steps 3 and 4 are for stochastic models and are the part no other tool does.
Background on what we have today
The implementation is in
pybnf/algorithms/optimizers/scatter_search.py. It keeps a small set of good parameter sets called the reference set. Each iteration it combines every ordered pair of reference members to produce candidate parameter sets, simulates all of them, and replaces a reference member when one of its children scores better. A member that fails to improve for several iterations is retired to an archive and replaced with a fresh random point.Against Glover's original description we have the diverse starting population, the systematic combination of every pair, and the combination formula itself, which comes from Egea and colleagues. Two pieces are missing.
Step 1: add the improvement method
Glover's template applies a local search to the candidates produced by combination. Egea's version does the same. Ours does no local search at all and relies purely on recombination. This is the largest gap and the most likely source of better fits.
We already have the part needed.
SimplexRunnerinpybnf/algorithms/optimizers/simplex.pyis a self-contained, picklable state machine that runs one Nelder-Mead simplex search, works in parameter set space, and knows nothing about the trajectory, backup, or the job scheduler.PowellRunnerinpybnf/algorithms/optimizers/powell.pyis the same shape. Both were built for the concurrent multi-start work and can be driven from inside scatter search.Local search is expensive, so it should not run on every candidate. Egea's version uses filters to decide when a candidate is promising enough to be worth refining. We need something similar. On a cluster the refinements can run at the same time as each other, so the cost is wall clock only when we are short of processors.
Step 2: choose the diverse half of the reference set by distance
round_1_initat line 146 builds the first reference set by taking the best half by score and then filling the rest withrng.choiceover the remainder. Glover's point is that the second half should be the most diverse members, chosen by their distance from what is already in the set. Picking at random gives diversity only on average, and in a problem with many parameters that is much weaker than choosing for it deliberately.This is a small change to one method and should be measurable on its own.
Step 3: make the reference set aware of noise
When a model is stochastic, running it twice with the same parameters gives two different objective values. Scatter search currently stores one simulation's result as a member's score and treats it as fact.
The useful observation is that scatter search never needs the objective values to be precise. It only needs the ordering to be right. The step size in the combination formula comes from the gap in rank between the two members being combined:
Sorting the reference set, deciding whether a child beat its parent, and counting whether a member is stuck are all ranking decisions too.
So the reference set should hold an estimate and a measure of uncertainty for each member rather than a single number, and extra simulations should be spent only where the ordering is genuinely in doubt. Two members whose uncertainty ranges overlap, and whose rank gap is driving a combination, are worth more simulations. Members that are clearly separated are worth none. How many repeats a parameter set receives then follows from how contested its position is, rather than from a setting the user has to choose.
This also fixes a real fault. The stuck counter at line 197 does not know about noise. A member that got a lucky simulation cannot be beaten by honest children, so its counter climbs until it is retired into the
local_minsarchive as a local minimum that it never was. The archive of best results quietly fills with artifacts of the noise.Related work is in #659, which fixes the same kind of error in the final answer at the end of a run. That issue is smaller, it stands alone, and it will tell us how large this problem actually is on real fits before we commit to this step.
Step 4: stop waiting for the whole iteration, and give idle processors something useful to do
ScatterSearchsetswaits_for_full_generation = True, so every iteration waits for its slowest simulation. This is worst for exactly the models we care about here, because the running time of a stochastic simulation depends on the trajectory it happens to take, so the spread in running times is wide.Once step 3 exists, the fix for the idle processors is already in hand. When a processor becomes free and no new candidate is ready, give it another repeat of a parameter set whose rank is still uncertain. The idle processor problem and the noise problem then share one answer.
How this differs from saCeSS
Penas and colleagues published a parallel, cooperative, self-adapting scatter search in 2017, and our source already notes that we implemented the plain version and not theirs. It is worth being clear about why this epic is not a copy of it.
Their method is built for large deterministic models, where the objective is a fixed number. A processor can either evaluate a new candidate or run another independent copy of the search. With a stochastic model there is a third choice, which is to reduce the uncertainty about a candidate already seen, and their design has no way to express it because it has no notion of noise.
They also run many independent copies of the search, each doing its combinations one after another. We run a single search and evaluate all of its combinations at once. Both are reasonable uses of a large machine, and ours is already built.
Also worth doing
The size of the reference set decides how many simulations run at once, so the right value is close to the square root of the number of processors available. We already read the number of connected workers in
Algorithm._report_parallelism. We could set a sensible reference set size automatically instead of asking users to work outpopulation_size * (population_size - 1)against their allocation by hand, which is whatdocs/cluster.rsttells them to do today.