Apply pykrige to interpolate a high-dimensional data

I am going to realize a black optimization problem, where one step is applying Kriging interpolation method to construct a approximate model. I have referred to the PyKrige Documentation and do not find any class supporting interpolation of high-dimensional(>3) dataset.

The following describes a 4-dim example:

from pykrige.uk import UniversalKriging
import numpy as np

data = np.array([[0.3, 0.6, 1.1, 1.2, 0.47],
                [1.9, 1.3, 0.4, 0.6, 0.56],
                [1.1, 1.0, 2.9, 3.2, 0.74],
                [3.3, 3.1, 4.4, 4.5, 1.47],
                [4.7, 4.3, 3.8, 3.9, 1.74]])

UK = UniversalKriging(x=data[:, 0:2], y=data[:, 2:4], z=data[:, 4], variogram_model='linear', drift_terms=['regional_linear'])

However, it reports an error "operands could not be broadcast together with shapes (2,10) (1,2) (2,10) ".

So, which class should I use to realize high-dimensional dataset?

You could use GSTools, the more general sister package of PyKrige:

import gstools as gs
import numpy as np

data = np.array(
    [
        [0.3, 0.6, 1.1, 1.2, 0.47],
        [1.9, 1.3, 0.4, 0.6, 0.56],
        [1.1, 1.0, 2.9, 3.2, 0.74],
        [3.3, 3.1, 4.4, 4.5, 1.47],
        [4.7, 4.3, 3.8, 3.9, 1.74],
    ]
)
# choose your model wisely (linear is actually not valid in 4D)
model = gs.Linear(dim=4, var=1, len_scale=10)

x0 = data[:, 0]
x1 = data[:, 1]
x2 = data[:, 2]
x3 = data[:, 3]
c = data[:, 4]

# unstructured output grid at 10000 random locations
seed = gs.random.MasterRNG(19970221)
rng = np.random.RandomState(seed())
gx0 = rng.uniform(0, 5, size=10000)
gx1 = rng.uniform(0, 5, size=10000)
gx2 = rng.uniform(0, 5, size=10000)
gx3 = rng.uniform(0, 5, size=10000)
grid = (gx0, gx1, gx2, gx3)

ok = gs.krige.Ordinary(model, cond_pos=(x0, x1, x2, x3), cond_val=c)
field, krige_var = ok(grid, mesh_type="unstructured")

ok.plot()