The accurate prediction of the compression index (Cc) is crucial for understanding the settlement behavior of clayey soils, which is a key factor in geotechnical design. Traditional empirical models, while widely used, often fail to generalize across diverse soil conditions due to their reliance on simplified assumptions and regional dependencies. This study proposed a novel hybrid method combining Genetic Programming (GP) and XGBoost methods. A large database (including 385 datasets) of geotechnical properties, including the liquid limit (LL), the plasticity index (PI), the initial void ratio (e0), and the water content (w), was used. The hybrid GP-XGBoost model achieved remarkable predictive performance, with an R2 of 0.966 and 0.927 and mean squared error (MSE) values of 0.001 and 0.001 for training and testing datasets, respectively. The mean absolute error (MAE) was also exceptionally low at 0.030 for training and 0.028 for testing datasets. Comparative analysis showed that the hybrid model outperformed the standalone GP (R2 = 0.934, MSE = 0.003) and XGBoost (R2 = 0.939, MSE = 0.002) models, as well as traditional empirical methods such as Terzaghi and Peck (R2 = 0.149, MSE = 0.090). Key findings highlighted that the initial void ratio and water content are the most influential predictors of Cc, with feature importance scores of 0.55 and 0.27, respectively. The novelty of the proposed method lies in its ability to combine the interpretability of GP with the computational efficiency of XGBoost and results in a robust and adaptable predictive tool. This hybrid approach has the potential to advance geotechnical engineering practices by providing accurate and interpretable models for diverse soil profiles and complex site conditions.