Efficient and accurate simulations of wave propagation are central to applications in seismology. In practice, heterogeneities arise from the presence of different types of rock in the subsurface. Additionally, simulations over long time periods require high order approximation to avoid numerical dispersion and dissipation effects. The weight-adjusted discontinuous Galerkin (WADG) method delivers high order accuracy for arbitrary heterogeneous media. However, the cost of WADG grows rapidly with the order of approximation. To reduce the computational complexity of high order methods, we propose a Bernstein-Bézier WADG method, which takes advantage of the sparse structure of matrices under the Bernstein-Bézier basis. Our method reduces the computational complexity from O(N^6) to O(N^4) in three dimensions and is highly parallelizable to implement on Graphics Processing Units (GPUs).