We address the problem of efficient phase diagram sampling by adopting active learning techniques from machine learning and achieve drastic reduction in the sample size (number of sampled state points) needed to establish the phase boundary up to a given precision. Although advanced sampling techniques such as Gibbs ensemble and Gibbs-Duhem integration can sample phase equilibria efficiently, they may fail to generalize to many nonequilbrium systems. This forces researchers to resort to grid search simulations when studying many important active matter systems. Grid search suffers from low efficiency by sampling predetermined state points that provide little information about the phase boundaries. We propose an active learning framework to overcome this deficiency by adaptively choosing the next most informative state point(s) every round. This is done by interpolating the sampled state points' phases by Gaussian process regression and then using an acquisition function to quantify the informativeness of possible next state points. We generalize our approach with asynchronous sampling techniques to better utilize parallel computing resources and extend the algorithm to incorporate uncertainty information from multiple replicas. We demonstrate the usefulness of our approach in two example systems in soft matter physics: a phase-separated steady-state of a mixture of self-propelled and passive Brownian colloids, and an equilibrium phase boundary between two quasicrystals. We achieve 5 times sample efficiency improvement in the case of the former example while the phase diagram in the latter example has not been studied before with comparable precision. Our active learning enhanced phase diagram sampling method greatly accelerates research and opens up opportunities for extra-large-scale exploration of a wide range of phase diagrams by simulation or experiment.