  We study a family of regularized score-based estimators for learning the
structure of a directed acyclic graph (DAG) for a multivariate normal
distribution from high-dimensional data with $p\gg n$. Our main results
establish support recovery guarantees and deviation bounds for a family of
penalized least-squares estimators under concave regularization without
assuming prior knowledge of a variable ordering. These results apply to a
variety of practical situations that allow for arbitrary nondegenerate
covariance structures as well as many popular regularizers including the MCP,
SCAD, $\ell_{0}$ and $\ell_{1}$. The proof relies on interpreting a DAG as a
recursive linear structural equation model, which reduces the estimation
problem to a series of neighbourhood regressions. We provide a novel
statistical analysis of these neighbourhood problems, establishing uniform
control over the superexponential family of neighbourhoods associated with a
Gaussian distribution. We then apply these results to study the statistical
properties of score-based DAG estimators, learning causal DAGs, and inferring
conditional independence relations via graphical models. Our results
yield---for the first time---finite-sample guarantees for structure learning of
Gaussian DAGs in high-dimensions via score-based estimation.
