GPflow 2.0 中的自定义 Haversine Matern52 内核

Question

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

当前编写的class不会改变Matern52 class；使用此 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 和 develop 分支中，对于 Stationary 的子 classes，K 调用 self.scaled_squared_euclid_dist -- 但是您在 Haversine 版本中定义的方法称为 scaled_squared_dist，所以这是一个 new 方法，您实际上并没有从 Matern52 内核 class 覆盖它的基础 class 方法！（也许 gpflow.kernels.stationaries.Stationary 中的方法最好称为 scaled_squared_dist。）

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

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

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 距离内起作用 space 因此需要适合它 - 取决于您的数据集。在构造基本内核时将其作为 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)

GPflow 2.0 中的自定义 Haversine Matern52 内核

Custom Haversine Matern52 kernel in GPflow 2.0

python

tensorflow

tensorflow-probability

gpflow

tensorflow2.0