GPflow 2.0 中的自定义 Haversine Matern52 内核

Question

使用 GPflow 2.0，我想用 Haversine 而不是欧几里德距离实现自定义的 Matern 5/2 内核。我在gpflow.kernels.Matern52 类之上创建了一个自定义类，其中包含一个scaled_squared_dist 函数来覆盖从 Stationary 类继承的scaled_squared_euclid_dist。

当前编写的类不会改变 Matern52 类；使用 HaversineKernel_Matern52 内核的 GP 回归与使用 Matern52 内核的 GP 回归完全一样。

import gpflow
from gpflow.utilities.ops import square_distance

class HaversineKernel_Matern52(gpflow.kernels.Matern52):
    """
    Isotropic Matern52 Kernel with Haversine distance instead of euclidean distance.
    Assumes 2-dimensional inputs, with columns [latitude, longitude] in degrees.
    """

    def __init__(self, lengthscale=1.0, variance=1.0, active_dims=None, ard=None):
        super().__init__(active_dims=active_dims, variance=variance, 
                         lengthscale=lengthscale, ard=ard)

    def haversine_dist(self, X, X2):
        pi = np.pi / 180
        f = tf.expand_dims(X * pi, -2)  # ... x N x 1 x D
        f2 = tf.expand_dims(X2 * pi, -3)  # ... x 1 x M x D
        d = tf.sin((f - f2) / 2) ** 2
        lat1, lat2 = tf.expand_dims(X[:, 0] * pi, -1), \
                    tf.expand_dims(X2[:, 0] * pi, -2)
        cos_prod = tf.cos(lat2) * tf.cos(lat1)
        a = d[:,:,0] + cos_prod * d[:,:,1]
        c = tf.asin(tf.sqrt(a)) * 6371 * 2
        return c

    def scaled_squared_dist(self, X, X2=None):
        """
        Returns ||(X - X2ᵀ) / ℓ||² i.e. squared L2-norm.
        """

        X_scaled = X / self.lengthscale
        X2_scaled = X2 / self.lengthscale if X2 is not None else X2
        return square_distance(X_scaled, X2_scaled)

我需要改变什么才能使这个内核正确地重新计算Haversine距离？

这个问题建立在GPflow issue #969 之上。

谢谢！

Answer 1

GP 代码使用内核的K（和K_diag）方法。在 GPflow 2.0.0rc1 和开发分支中，对于 Stationary 的子类，K 调用 self.scaled_squared_euclid_dist -- 但是您在 Haversine 版本中定义的方法称为 scaled_squared_dist，所以这是一个 new 方法，您实际上并没有从Matern52 内核类覆盖它的基类方法！ （也许 gpflow.kernels.stationaries.Stationary 中的方法最好称为scaled_squared_dist。）

此外，您的scaled_squared_dist 只调用square_distance 而不是self.haversine_dist；假设后者返回一个距离，而不是它的平方，您还需要将它包装在tf.square() 中。 haversine_dist 方法似乎也没有考虑到lengthscale 参数。

如果您想尝试几个具有Haversine 距离的不同内核，一种更健壮/可重用的编码方式可能是编写一个将任何固定内核作为参数的包装类，并重新定义核矩阵方法：

def haversine_dist(X, X2):
    pi = np.pi / 180
    f = tf.expand_dims(X * pi, -2)  # ... x N x 1 x D
    f2 = tf.expand_dims(X2 * pi, -3)  # ... x 1 x M x D
    d = tf.sin((f - f2) / 2) ** 2
    lat1, lat2 = tf.expand_dims(X[:, 0] * pi, -1), \
                tf.expand_dims(X2[:, 0] * pi, -2)
    cos_prod = tf.cos(lat2) * tf.cos(lat1)
    a = d[:,:,0] + cos_prod * d[:,:,1]
    c = tf.asin(tf.sqrt(a)) * 6371 * 2
    return c


class HaversineDistance(gpflow.kernels.Stationary):
    def __init__(self, base: gpflow.kernels.Stationary):
        self.base = base

    @property
    def active_dims(self):
        return self.base.active_dims

    @property
    def variance(self):
        return self.base.variance  # for K_diag to work

    def K(self, X, X2=None, presliced=False):
        if not presliced:
            X, X2 = self.slice(X, X2)
        if X2 is None:
            X2 = X

        # Note that lengthscale is in Haversine distance space:
        r = haversine_dist(X, X2) / self.base.lengthscale

        if hasattr(self.base, "K_r"):
            return self.base.K_r(r)
        else:
            return self.base.K_r2(tf.square(r))

我将您的haversine_dist 分解为它自己的函数（这将使编写测试更容易！）。这适用于任何固定内核，无论是SquaredExponential（定义K_r2）还是Matern52（定义K_r）。你可以简单地使用它：

k = HaversineDistance(gpflow.kernels.SquaredExponential())

或使用 Matern 3/2 作为基础内核：

k = HaversineDistance(gpflow.kernels.Matern32())

请注意，与往常一样，超参数初始化很重要，您需要记住，基本内核的长度尺度在Haversine 距离空间内起作用，因此需要适合它 - 取决于您的数据集。在构建基本内核时将其作为lengthscale 和variance 参数传入：

k = HaversineDistance(gpflow.kernels.Matern32(variance=5.0, lengthscale=0.1))

或者在基础内核的 gpflow.Parameter 属性上使用assign()：

k.base.variance.assign(5.0)
k.base.lengthscale.assign(0.1)