Dear UQLab developers,
with this post, I want to highlight what I believe to be a possible bug I encountered with the Kriging module. Specifically, the correlation function lengthscale “theta” returned in the output structure (Model.Kriging.theta) does not match the actual and correct value that is computed internally.
To support my claim, let me consider the simple test case considered in Section 2.3 of the user manual of the Kriging module. The code reported at the end of this post implements the simulation, with minimal and trivial modifications w.r.t. the code available in the manual. Please note that, differently from the manual’s example, I use a Gaussian correlation function instead of the default Matern 5/2 function. However, that the same problem is found with the Matern correlation as well.
The output structure returns a value of theta = 0.5276. My custom definition of the leave-one-out cross-validation function to estimate the lengthscale (not reported in this post for brevity) suggests instead a value of theta = 2.5877.
I believe that the values I computed are correct, and I provide two evidences to support my claim. The first is that the correlation matrix R of the training data in the output structure does not match the one that is generated by evaluating the correlation function with the returned lengthscale, as is verified by running the code I provide.
The second proof is obtained by creating a manual Kriging predictor (prediction-only mode) with the UQLab-returned parameters and comparing the predictions. Note that the predictions of the two models match only if “my_theta” is used in place of the UQLab-returned theta.
Note that the same problem is found with the default Matern 5/2 correlation function (in that case, returned and correct values of theta are 2.9060 and 14.2361, respectively).
This are the UQLab / MATLAB / OS versions I’m using:
UQLab version: UQLab_Rel2.0.0
UQLab modules in question: Kriging
MATLAB version: R2023a
Operating system: Windows 10 Pro
Code:
% Start the UQLab framework
uqlab
rng(100,'twister') % For reproducible results
% Define function
fun = @(X) X.*sin(X);
% Plot function
X = transpose(0:0.01:14);
Y = fun(X);
figure(1)
plot(X,Y,'LineWidth',1.5); hold on;
% Create experimental design
X_train = transpose(0:2:14);
Y_train = fun(X_train);
% Plot training samples
plot(X_train,Y_train,'ko','MarkerFaceColor','k');
% Define a Kriging metamodel
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'Kriging';
MetaOpts.ExpDesign.Sampling = 'User';
MetaOpts.ExpDesign.X = X_train;
MetaOpts.ExpDesign.Y = Y_train;
MetaOpts.Corr.Family = 'gaussian';
% Create the Kriging metamodel
myKriging = uq_createModel(MetaOpts);
% Evaluate the Kriging metamodel
Y_pred = uq_evalModel(myKriging,X);
figure(1)
plot(X,Y_pred,'--','LineWidth',2.0);
% Correlation matrix computed with UQLab-returned lengthscale
disp('Correlation matrix computed with UQLab-returned lengthscale:');
R = uq_eval_Kernel(X_train,X_train,myKriging.Kriging.theta,myKriging.Internal.Kriging.GP.Corr)
% Compare with UQLab correlation matrix
disp('Compare with UQLab correlation matrix:');
myKriging.Internal.Kriging.GP.R
% Correlation matrix computed with my lengthscale
disp('Correlation matrix computed with my lengthscale:');
my_theta = 2.5877;
my_R = uq_eval_Kernel(X_train,X_train,my_theta,myKriging.Internal.Kriging.GP.Corr)
% Comparison in terms of relative error norm
disp('Comparison in terms of relative error norm:');
norm(R-myKriging.Internal.Kriging.GP.R)/norm(myKriging.Internal.Kriging.GP.R)
norm(my_R-myKriging.Internal.Kriging.GP.R)/norm(myKriging.Internal.Kriging.GP.R)
% Define a Kriging predictor-only model with UQLab-returned lengthscale
clear MetaOpts;
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'Kriging';
MetaOpts.ExpDesign.Sampling = 'User';
MetaOpts.ExpDesign.X = X_train;
MetaOpts.ExpDesign.Y = Y_train;
% Select ordinary Kriging
MetaOpts.Kriging.Trend.Type = 'ordinary';
% Set the values of beta , sigmaˆ2 and theta
MetaOpts.Kriging.beta = myKriging.Kriging.beta;
MetaOpts.Kriging.sigmaSQ = myKriging.Kriging.sigmaSQ;
MetaOpts.Kriging.theta = myKriging.Kriging.theta; % uses UQLab-returned lengthscale
%MetaOpts.Kriging.theta = my_theta; % uses lengthscale computed by me
% Select ellipsoidal, Gaussian correlation function
MetaOpts.Kriging.Corr.Type = 'ellipsoidal';
MetaOpts.Kriging.Corr.Family = 'gaussian';
% Create the metamodel
myCustomKriging = uq_createModel(MetaOpts);
% Evaluate the metamodel on some new points
Y_pred = uq_evalModel(X);
figure(1)
plot(X,Y_pred,':','LineWidth',2.5);