The Hawkes process has been widely applied to modeling self-exciting events, including neuron spikes, earthquakes and tweets. To avoid designing parametric kernel functions and to be able to quantify the prediction confidence, non-parametric Bayesian Hawkes processes have been proposed. However the inference of such models suffers from unscalability or slow convergence. In this paper, we first propose a new non-parametric Bayesian Hawkes process whose triggering kernel is modeled as a squared sparse Gaussian process. Second, we present the variational inference scheme for the model optimization, which has the advantage of linear time complexity by leveraging the stationarity of the triggering kernel. Third, we contribute a tighter lower bound than the evidence lower bound of the marginal likelihood for the model selection. Finally, we exploit synthetic data and large-scale social media data to validate the efficiency of our method and the practical utility of our approximate marginal likelihood. We show that our approach outperforms state-of-the-art non-parametric Bayesian and non-Bayesian methods.