For complex nonlinear systems, it is challenging to design algorithms that are fast, scalable, and give an accurate approximation of the stability region. This paper proposes a sampling-based approach to address these challenges. By extending the parametrization of quadratic Lyapunov functions with the system dynamics and formulating an $\ell_1$ optimization to maximize the invariant set over a grid of the state space, we arrive at a computationally efficient algorithm that estimates the domain of attraction (DOA) of nonlinear systems accurately by using only linear programming. The scalability of the Lyapunov function synthesis is further improved by combining the algorithm with ADMM-based parallelization. To resolve the inherent approximative nature of grid-based techniques, a small-scale nonlinear optimization is proposed. The performance of the algorithm is evaluated and compared to state-of-the-art solutions on several numerical examples.