  This article focuses on numerical issues in maximum likelihood parameter
estimation for Gaussian process regression (GPR). This article investigates the
origin of the numerical issues and provides simple but effective improvement
strategies. This work targets a basic problem but a host of studies,
particularly in the literature of Bayesian optimization, rely on off-the-shelf
GPR implementations. For the conclusions of these studies to be reliable and
reproducible, robust GPR implementations are critical.
