For machine learning of interatomic potentials a scalable sparse Gaussian process regression formalism is introduced with a data-efficient on-the-fly adaptive sampling algorithm. With this approach, the computational cost is effectively reduced to those of the Bayesian linear regression methods whilst maintaining the appealing characteristics of the exact Gaussian process regression. As a showcase, experimental melting and glass-crystallization temperatures are reproduced for Li7P3S11, Li diffusivity is simulated, and an unchartered phase is revealed with much lower Li diffusivity which should be circumvented.