Directed acyclic graphs (DAGs) are widely used to describe directional relations among interacting units. Directional relations are estimated by reconstructing a DAG's structure, which is a great challenge when the total ordering of a DAG is unknown. In such a situation, existing methods such as the neighborhood and search-and-score methods suffer greatly, as the overall estimation error accumulates super-exponentially in the number of nodes, especially when a local/sequential approach enumerates edge directions by testing or optimizing a criterion locally. In other words, a local method may break down even for moderately sized graphs. In this thesis, we propose a novel approach to simultaneously identify all estimable directed edges as well as model parameters jointly. This approach uses constrained maximum likelihood with nonconvex constraints reinforcing acyclicity. Computationally, we develop a novel reduction method that constructs a set of active constraints (cubic in the number of nodes) from the super-exponentially many constraints. This, coupled with an alternating direction method of multipliers and a difference convex method, permits efficient computation for large graph learning. Theoretically, we show that the proposed method consistently reconstructs identifiable directions of the true graph, under a degree of reconstructability assumption. This goes beyond the strong faithfulness assumption, commonly used in the literature. Moreover, the method recovers the optimal performance of the oracle estimator in terms of parameter estimation. Numerically, the method compares favorably against its competitors.