gpflow中具有减法核的Cholesky分解问题

Question

我正在尝试以下内核：

class AperiodicMatern12(gpflow.kernels.Kernel):
  def __init__(self, input_dim = None, period = 1.0, variance = 1.0, 
      lengthscales = 1.0, active_dims = None, name = None):
    super().init(input_dim, active_dims = active_dims)
    k0 = gpflow.kernels.Matern12(input_dim, active_dims = active_dims)
    k1 = gpflow.kernels.Matern12(input_dim, active_dims = active_dims)
    k = gpflow.kernels.Periodic(base = k1, period = period)
    self.base = k0
    self.per = k

  @gpflow.params_as_tensors
  def K(self, X, X2 = None):
    res = self.base.K(X, X2) - self.per.K(X, X2)
    return res

  def Kdiag(self, X):
    return np.diag(self.K(X))

当我运行它时：

lik = gpflow.likelihoods.Gaussian()
k = AperiodicMatern12(1, active_dims = [0])
m = gpflow.models.GPR(X, Y, kern = k)
gpflow.train.ScipyOptimizer().minimize(m)

我得到：

InvalidArgumentError (回溯见上文)：Cholesky分解没有成功。输入可能无效。[节点GPR-0b2840db-15/似然_1/Cholesky(定义为/Users/mjg/anaconda3/lib/python3.6/site-packages/gpflow/models/gpr.py:72) )

当我尝试使用定义为加法的内核时：

@gpflow.params_as_tensors
def K(self, X, X2 = None):
  res = self.base.K(X, X2) + self.per.K(X, X2)
  return res

一切都很好。X和Y在这两种情况下都是标准化的。那么减法内核有什么问题呢？

Answer 1

Cholesky分解只适用于正定矩阵.所以Cholesky分解操作会失败，如果输入参数不是正定的，或者--换句话说--当输入矩阵有任何负的特征值时。当这种情况发生时，实际计算核矩阵(例如np.linalg.eigvals(tf.Session().run(k.K(X))))的特征值通常是有帮助的。

核矩阵的特征值可能是负的有两个不同的情况:由于有限的数值精度问题，最常见的原因是“几乎正”负特征值--这就是为什么在计算GPflow中的Cholesky分解之前，我们通常添加一个抖动矩阵(对角线上有10^(-6)的对角矩阵)。第二个原因是核矩阵的计算存在误差。这里的情况如下：

X = np.linspace(0, 10, 11)[:,None]  # dummy 1D data matrix
session = gpflow.get_default_session()
K = session.run(k.K(X))
np.linalg.eigvals(K)

输出

array([-9.00767547,  1.78835489,  1.51465032,  1.18585391,  0.97478308,
        0.79231475,  0.67625953,  0.5866432 ,  0.52938737,  0.4691057 ,
        0.49032271])

虽然两个核(正定矩阵)的和总是给出一个有效的核(另一个正定矩阵)，但这种性质在一般情况下不成立。

更新:在非常特殊的情况下，两个内核的差异可能是另一个有效的正定内核。这就是Durrande等人关于周期和非周期母类核的构造的情况。它们的构造是特定于母体家族的内核的，并且与基于在gpflow.kernels.Periodic中实现的任何固定内核的周期内核的一般实现不同。

GPflow不提供Durrande等人。虽然在GPflow上实现它们很容易，但是GPflow项目很乐意考虑将这些内核添加到GPflow代码库中的拉请求。