  It is well known in the literature that the problem of learning the structure
of Bayesian networks is very hard to tackle: its computational complexity is
super-exponential in the number of nodes in the worst case and polynomial in
most real-world scenarios.
  Efficient implementations of score-based structure learning benefit from past
and current research in optimisation theory, which can be adapted to the task
by using the network score as the objective function to maximise. This is not
true for approaches based on conditional independence tests, called
constraint-based learning algorithms. The only optimisation in widespread use,
backtracking, leverages the symmetries implied by the definitions of
neighbourhood and Markov blanket.
  In this paper we illustrate how backtracking is implemented in recent
versions of the bnlearn R package, and how it degrades the stability of
Bayesian network structure learning for little gain in terms of speed. As an
alternative, we describe a software architecture and framework that can be used
to parallelise constraint-based structure learning algorithms (also implemented
in bnlearn) and we demonstrate its performance using four reference networks
and two real-world data sets from genetics and systems biology. We show that on
modern multi-core or multiprocessor hardware parallel implementations are
preferable over backtracking, which was developed when single-processor
machines were the norm.
