Efficient monotonic Gaussian processes via Randomize-then-Optimize

• New optimization-based virtual point method RLRTO for monotonic Gaussian Processes • Open-source Python implementation of RLRTO and NUTS methods using the CUQIpy library • Extensive numerical experiments demonstrate superior efficiency of RLRTO over MCMC Gaussian processes (GPs) are widely used as surrogate models for complicated functions in scientific and engineering applications. In many cases, prior knowledge about the function to be approximated, such as monotonicity, is available and can be leveraged to improve model fidelity. Incorporating such constraints into GP models enhances predictive accuracy and reduces uncertainty, but remains a computationally challenging task for high-dimensional problems. In this work, we present a novel virtual point-based framework for building constrained GP models under monotonicity constraints, based on Regularized Linear Randomize-then-Optimize (RLRTO), which enables efficient sampling from a constrained posterior distribution by means of solving randomized optimization problems. We compare the proposed method with commonly used MCMC methods like Gibbs sampling and the No U-Turn Sampler (NUTS). A Python implementation of these methods is provided and can be easily applied to a wide range of problems. This implementation is then used to validate the proposed method on approximating a range of synthetic functions, demonstrating comparable predictive performance between all considered methods and significant improvements in computational efficiency with the proposed RLRTO method. The framework is further applied to construct surrogate models for systems of differential equations.