Shear constitutive models of rock discontinuities have been viewed as an effective stability evaluation tool in the rock mass engineering application area. This paper proposes a new interpretable shear strength criterion for rock joints based on multivariate adaptive regression splines (MARS) algorithm. Sensitivity analyses are then performed on the developed shear strength criterion. As the second purpose of this research, the potential competence of five surrogate soft computing (SC) methods including Gaussian process (GP), alternating model tree (AMT), Cubist, radial basis function (RBF) networks, and elastic net (EN) paradigms for fast predicting the shear strength of rock discontinuities is also comparatively evaluated along with the MARS model. These approaches formulate nonlinear relations between input and output variables. The proposed methodologies consider eight input factors: sampling interval (1), maximum contact area ratio (A(0)), distribution parameter (C), maximum apparent dip angle (theta(max)*), basic friction angle (phi(b)), tensile strength (sigma(t)), uniaxial compressive strength (sigma(c)), and normal stress (sigma(n)) to assess peak shear strength (tau(p)) of rock joints. A dataset collected from the literature conducted on the direct shear test is used for training and evaluating proposed methods. Ten-fold cross-validation is utilized to enhance the robustness and generalization of the developed SC-based shear strength surrogate models. Statistical indices and comparative analyses with conventional criteria indicate that the proposed SC regressors can deliver good agreements with the measured data in terms of performance accuracy and outperform or achieve comparable performances to the conventional models. However, GP, AMT, and RBF models have a good prediction performance than MARS and EN models.