1

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

当前编写的类不会改变 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 问题 #969之上。

谢谢！

4

1 回答 1

4

GP 代码使用内核的K(and K_diag) 方法。在 GPflow 2.0.0rc1 和开发分支中，对于 的子类Stationary，K调用self.scaled_squared_euclid_dist- 但是您在 Haversine 版本中定义的方法被调用scaled_squared_dist，所以这是一个新方法，您实际上并没有从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 距离空间内起作用，因此需要适合它 - 取决于您的数据集。在构建基本内核时将其作为lengthscaleand参数传入：variance

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

或者通过使用assign()基本内核的 gpflow.Parameter 属性：

k.base.variance.assign(5.0)
k.base.lengthscale.assign(0.1)
于 2019-12-28T22:36:42.977 回答