Efficient evaluation of slope stability is a frontier in geo‐disaster prevention fields. While many slope stability evaluation methods, ranging from deterministic to probabilistic, have been proposed, reliability methods are particularly advantageous as they can account for uncertainties during slope stability evaluation. To address the problem of the low efficiency of direct Monte Carlo simulations and overcome the defects of the traditional response surface method, in this work, a novel probabilistic procedure that integrates a Gaussian process regression‐based surrogate model and the random limit equilibrium method for slope stability evaluation while accounting for the spatial variability of soil strength is proposed. The novel Gaussian process regression‐based surrogate model is efficiently used in the Monte Carlo simulation to reduce the number of calls for direct stability analysis of a slope with spatially varied soil strength. To verify the accuracy and efficiency of the proposed procedure, it was applied to three case studies: the first one is a slope in saturated clay under undrained conditions; the second one is a slope for which the friction angle and cohesion values are cross‐correlated; and the third one is a real slope with multiple soil layers. The results obtained from the comparisons with other methods confirmed the precision and feasibility of the proposed procedure.