We develop a penalized likelihood approach to estimating the structure of a Gaussian Bayesian network, given by a directed acyclic graph, from observational data under a concave penalty. The framework introduced here does not rely on faithfulness or knowledge of the ordering of the variables and favours sparsity over complexity in estimating the underlying graph. Asymptotic theory for the estimator is provided in the finite-dimensional case, and a fast numerical scheme is offered that is capable of estimating the structure of graphs with thousands of nodes. By reparametrizing the usual normal log-likelihood, we obtain a convex loss function which accelerates computation of the proposed estimator. Our algorithm also takes advantage of sparsity and acyclicity by using coordinate descent, a computational approach which has recently become quite popular. Finally, we compare our method with the well-known PC algorithm, and show that our method is faster in general and does a significantly better job of handling small samples and very sparse networks. Our focus is on the Gaussian linear model, however, the framework introduced here can also be extended to non-Gaussian and non-linear designs, which is an attractive prospect for future applications.