Noiseless Gausian Proces regression not crossing the training points

I'm trying to perform a Gaussian Process Regression on the following series of data points (X and Y):

X         Y
0.00001  -6.8200000000000003
0.00014  -6.8073846153846160
0.00042  -6.6818333333333326
0.00200  0.5002400000000000
0.00545  -0.0163174311926606
0.00545  -0.0163174311926606
0.00545  -0.0163174311926606
0.00552  -0.0317916666666667
0.00604 -0.1288741721854305
0.00756  -0.2954391534391534
0.00987  -0.3925845997973657
0.02201  -0.2940145454545455
0.02618  -0.2296286259541984
0.04000  -0.1075512500000000
0.05737  -0.0801517421602787
0.05761  -0.0801819444444445
0.05762  -0.0801947916666667
0.05767  -0.0801204506065858
0.05771  -0.0801724436741768
0.05772  -0.0801856152512998
0.05773  -0.0801980935875217
0.05782  -0.0801769896193772
0.05783  -0.0801899653979239
0.05788  -0.0801169257340242
0.05789  -0.0801297063903282
0.05804  -0.0801881034482759
0.05948  -0.0800737815126051
0.07886  -0.0843219264892269
0.10000  -0.0866689000000000

I'm using GPy to optimize the hyperparameters while fixing the noise variance to 0:

GP_regression.           |                value  |  constraints  |  priors
rbf.variance             |   1.9537980554282994  |      +ve      |
rbf.lengthscale          |  0.47928874454359643  |      +ve      |
Gaussian_noise.variance  |                  0.0  |   +ve fixed   |

Then, I compute the values of Y* at each test point X* using the optimized lengthscale and variance from GPy. I normalize the data when optimizing the hyperparamters, so I normalize X* at the beginning of the function and undo the normalization at the end. I perform the inversion of the kernel matrix using Cholesky decomposition and add a very small jitter to handle matrices with linearly dependent columns:

L = cholesky(K + np.eye(len(X)) * 1e-11, lower=True)  # Add jitter for stability

# Solve for alpha (K * alpha = Y)
alpha_imag= cho_solve((L, True), Y)

def gp_surrogate_model_imag(x_star):

    x_star = np.atleast_1d(x_star)
    x_star = x_star.reshape(-1,1)
    # Normalize input
    x_star_norm = (x_star - X_mean_old) / X_std_old
    # Compute covariance vector k(X^*, X)
    k_star = kernel.K(x_star_norm, X)

    Y_surrogate = k_star @ alpha_imag

    # Convert back to original scale (denormalization)
    Y_surrogate = Y_surrogate * Y_std_old + Y_mean_old

    return Y_surrogate.item()

Since the noise variance is 0, the points of the surrogate model should coincide exactly with the training points, right?

However, when I compute the differences between Y and Y* at the training points, I obtain differences on the order of 10^-5.

Is this an issue of my implementation, or are such small differences expected even when the noise variance is fixed to 0 in GPy?