In order to accurately and efficiently forecast the cutting force of complex surface multi-axis machining, a cutting force model based on Gaussian process regression (GPR) algorithm is developed in this paper. Feature parameters of tool–workpiece engagement, such as tool-axis inclination angles and cutting width, which serve as input characteristic parameters of GPR model for cutting force prediction in complex surface machining, are extracted based on cutter location file (CLS). The training set of the GPR model are obtained using the mechanical force model where the tool–workpiece engagement is calculated by Boolean operations. Cutting force simulation software for complex surface machining is developed and the efficiency of the proposed GPR model is verified by comparing with the traditional force prediction model which adopts the Boolean operations to calculate the tool–workpiece engagement. The error of cutting force prediction is less than 10% and the evaluation coefficient of prediction results is maintained above 0.98. An impeller runner machining experiment was designed to verify the accuracy of the GPR model proposed in this paper for predicting cutting force in machining complex curved surface. In force prediction based on the same CLS file for an impeller passage processing, the method using Boolean operation takes 161 s, while the time elapsed of the proposed model is only 1.63 s. The results indicate that the proposed model is efficient and accurate for cutting force prediction in complex surface machining.