  Learning a Bayesian network (BN) from data can be useful for decision-making
or discovering causal relationships. However, traditional methods often fail in
modern applications, which exhibit a larger number of observed variables than
data points. The resulting uncertainty about the underlying network as well as
the desire to incorporate prior information recommend a Bayesian approach to
learning the BN, but the highly combinatorial structure of BNs poses a striking
challenge for inference. The current state-of-the-art methods such as order
MCMC are faster than previous methods but prevent the use of many natural
structural priors and still have running time exponential in the maximum
indegree of the true directed acyclic graph (DAG) of the BN. We here propose an
alternative posterior approximation based on the observation that, if we
incorporate empirical conditional independence tests, we can focus on a
high-probability DAG associated with each order of the vertices. We show that
our method allows the desired flexibility in prior specification, removes
timing dependence on the maximum indegree and yields provably good posterior
approximations; in addition, we show that it achieves superior accuracy,
scalability, and sampler mixing on several datasets.
