Should we standardize the data while doing Gaussian process regression? I am performing Gaussian process regression (GPR) and optimizing over hyper-parameters. I am using minFunc to perform all optimizations. My question is should we (or rather, can we) standardize the data before giving it to the objective function? If we do standardize, then the hyper-parameters will be learned according to the standardized data. However, at test time, assuming we get samples one-by-one, it won't be possible to standardize each sample indepdently, right? (Unless, we use some standardizing factors from the training data). If it matters, all the elements in my data are between -1 to 1, however, some columns may have a very small mean and variance as compared to the other columns.
So my question is, should we normalize the data while doing GPR?
P.S. Actually, I observe some weird behavior if I don't standardize my data. For example, minFunc suddenly gives me step direction is illegal error. Some online reading led me to believe that there is either a problem in your gradient calculation or your data is not standardized. I am sure about my gradient function calculation, I have also check it with the DerivativeCheck option. So, that leaves the possibility of data not being standardized. 
 A: I agree with Alexey Zaytsev's answer (and the discussion of that answer) that standardising is good to do for various reasons also on the outputs. However, I just want to add an example of why it can be important to standardise the outputs.
An example of where normalizing/standardizing the output is important is for small and noisy output values. The following figure illustrates it on a very simple example. MATLAB code below figure.
The blue dots represent noisy samples from a sin() function. The orange/red circles represent the predicted Gaussian Process values.


*

*In the top subplot the amplitude is unity (with noise).

*In the second subplot the output values have been scaled with 1e-5 and we can see that the Gaussian Process model (with default settings) predicts a constant model. The default settings optimise on the noise parameter and have a fairly high lower bound.

*In the third subplot the noise parameter is set to zero and not optimised over. The model over-fits in this case.

*The fourth subplot shows the standardised outputs and a model fitted on that data.

*In the last subplot the outputs are scaled back from the standardisation operation (not scaled back to the original values) and also the predicted values are scaled. Note that the predicted values are scaled using the mean and standard deviation from the training data standardisation.



function importance_normgp()
% small outputs

% data
x = 0:0.01:1; x = x(:);
xp = linspace(0.1, 0.9, length(x)); xp = xp(:);

% noise free model
y = sin(2*pi*x) + 5e-1*randn(length(x), 1);

% train and predict gp model
mdl = fitrgp(x, y);
yp = predict(mdl, xp);

figure
subplot(5, 1, 1)
plot(x, y, '.')
hold on
plot(xp, yp, 'o')
title('original problem')

%% make outputs small (below noise lower bound)
ym = y/1e5;

% train and predict gp model
mdlm = fitrgp(x, ym);
ypm = predict(mdlm, xp(:));

subplot(5, 1, 2)
plot(x, ym, '.')
hold on
plot(xp, ypm, 'o')
title('small outputs')

%% outputs small and set sigma = 0

% train and predict gp model
mdlm1 = fitrgp(x, ym, 'Sigma', 1e-12, 'ConstantSigma', true, 'SigmaLowerBound', eps);
ypm1 = predict(mdlm1, xp(:));

subplot(5, 1, 3)
plot(x, ym, '.')
hold on
plot(xp, ypm1, 'o')
title('small outputs and sigma = 0')

%% normalise/standardise
nu = mean(ym);
sigma = std(ym);
yms = (ym - nu)/sigma;

% train and predict gp model
mdlms = fitrgp(x, yms);
ypms = predict(mdlms, xp(:));

subplot(5, 1, 4)
plot(x, yms, '.')
hold on
plot(xp, ypms, 'o')
title('standardised outputs')

% rescale
ypms2 = ypms*sigma + nu;

subplot(5, 1, 5)
plot(x, ym, '.')
hold on
plot(xp, ypms2, 'o')
title('scaled predictions')

legend('true model', 'prediction', 'Location', 'best')
end

A: Yes, it is desirable to standardize the data while learning Gaussian processes regression. There are a number of reasons:

*

*In common Gaussian processes regression model we suppose that output $y$ has zero mean, so we should standardize $y$ to match our assumption.

*For many covariance function we have scale parameters in covariance functions. So, we should standardize inputs to get better estimation of parameters of covariance functions.

*Gaussian processes regression is prone to numerical problems as we have to inverse ill-conditioned covariance matrix. To make this problem less severe, you should standardize your data.

Some packages do this job for you, for example GPR in sklearn has an option normalize for normalization of inputs, while not outputs; see this.
A: I want to add/comment to Alexey Zaytsev's answer:

*

*For the outputs $y$ - we don't need to assume 0 mean, it just simplifies the calculations when going to the conditional distribution of the new (predicted/test) data given train data. Also, we only need to de-mean, we don't have to scale.


*For the inputs $x$- the kernels invoke a measure of distance between the different sample points. E.g., RBF looks at $\Vert x-x^\prime\Vert$ - it doesn't make sense to look at the distance when there are different scales, so you need to normalize the input data first (unless they are on the same scale already).
