In complex and unknown processes, global models are fitted over the entire input domain but often tend to perform poorly whenever the response surface exhibits non-stationary behavior and varying smoothness. A common approach is to use local models, which requires partitioning the input domain into subdomains and training multiple models, thereby adding significant complexity. Recognizing this limitation, this study addresses the need for models that represent the input–output relationship consistently over the full domain while still adapting to local variations in the response. It introduces a novel machine learning approach: the Polynomial Chaos Expanded Gaussian Process (PCEGP), leveraging polynomial chaos expansion to calculate input-dependent hyperparameters of the Gaussian process (GP). This provides a mathematically interpretable approach that incorporates non-stationary covariance functions and heteroscedastic noise estimation to generate locally adapted models. The model performance is compared to different algorithms in benchmark tests for regression tasks. The results demonstrate low prediction errors of the PCEGP, highlighting model performance that is often competitive with or better than previous methods. A key advantage of the presented model is its interpretable hyperparameters along with training and prediction runtimes comparable to those of a standard GP.