Abstract: The accurate prediction and modeling of processing residual stress is the foundation for the prediction and control of residual stress deformation. To address the challenge of accurately predicting surface residual stress in multi-axis milling of Ti-6Al-4V titanium alloy, a hybrid approach integrating Gaussian process regression (GPR) and radial basis function (RBF) neural network is proposed. The Bayesian uncertainty quantification capability of GPR is exploited to optimize the initial weights and thresholds of the RBF network, thereby reducing prediction variance and enhancing generalization stability under small-sample conditions. Based on a multi-factor orthogonal experimental design, surface residual stress data corresponding to various combinations of process parameters within the designed parameter range were obtained using X-ray diffraction technique. The GPR algorithm was employed to collaboratively optimize the initial weights and thresholds of the RBF neural network, and the network was subsequently trained to establish a nonlinear mapping relationship between process parameters and residual stresses. The validation results show that the mean prediction errors in the σ_x and σ_y directions are 6.86% and 10.12%, respectively, with a single prediction time of only 3.82 s. Compared with other mainstream prediction algorithms, the proposed method achieves higher prediction accuracy and faster computational speed for surface residual stress.