0

我需要构建一个函数来给出高斯过程的后验协方差。这个想法是使用 GPytorch 训练 GP,然后获取学习的超参数,并将它们传递给我的内核函数。(由于几个原因,我不能直接使用 GPyTorch)。

现在的问题是我无法重现预测。这里是我写的代码。我整天都在努力,但我找不到问题所在。你知道我做错了什么吗?

        from gpytorch.mlls import ExactMarginalLogLikelihood 
        import numpy as np
        import gpytorch
        import torch
        train_x1 = torch.linspace(0, 0.95, 50) + 0.05 * torch.rand(50)
        train_y1 = torch.sin(train_x1 * (2 * np.pi)) + 0.2 * torch.randn_like(train_x1)

        n_datapoints = train_x1.shape[0]

        def kernel_rbf(x1, x2, c, l):
            # my RBF function
            if x1.shape is ():
                x1 = np.atleast_2d(x1)
            if x2.shape is ():
                x2 = np.atleast_2d(x2)
            return c * np.exp(- np.matmul((x1 - x2).T, (x1 - x2)) / (2 * l ** 2))

        class ExactGPModel(gpytorch.models.ExactGP):
            def __init__(self, train_x, train_y, likelihood):
                super().__init__(train_x, train_y, likelihood)

                lengthscale_prior = gpytorch.priors.GammaPrior(3.0, 6.0)
                outputscale_prior = gpytorch.priors.GammaPrior(2.0, 0.15)

                self.mean_module = gpytorch.means.ConstantMean()
                self.covar_module = gpytorch.kernels.ScaleKernel(
                    gpytorch.kernels.RBFKernel(lengthscale_prior=lengthscale_prior),
                    outputscale_prior=outputscale_prior)

            def forward(self, x):
                mean_x = self.mean_module(x)
                covar_x = self.covar_module(x)
                return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

        likelihood = gpytorch.likelihoods.GaussianLikelihood()
        model = ExactGPModel(train_x1, train_y1, likelihood)

        # Find optimal model hyperparameters
        model.train()
        likelihood.train()

        mll = ExactMarginalLogLikelihood(likelihood, model)

        # Use the Adam optimizer
        optimizer = torch.optim.Adam(model.parameters(), lr=0.1)  # Includes GaussianLikelihood parameters
        training_iterations = 50
        for i in range(training_iterations):
            optimizer.zero_grad()
            output = model(*model.train_inputs)
            loss = -mll(output, model.train_targets)
            loss.backward()
            print('Iter %d/%d - Loss: %.3f' % (i + 1, training_iterations, loss.item()))
            optimizer.step()

        # Get the learned hyperparameters 
        outputscale = model.covar_module.outputscale.item()
        lengthscale = model.covar_module.base_kernel.lengthscale.item()
        noise = likelihood.noise_covar.noise.item()

        train_x1 = train_x1.numpy()
        train_y1 = train_y1.numpy()

        # Get covariance train points
        K = np.zeros((n_datapoints, n_datapoints))
        for i in range(n_datapoints):
            for j in range(n_datapoints):
                K[i, j] = kernel_rbf(train_x1[i], train_x1[j], outputscale, lengthscale)

        # Add noise
        K += noise ** 2 * np.eye(n_datapoints)

        # Get covariance train-test points
        x_test = torch.rand(1, 1)
        Ks = np.zeros((n_datapoints, 1))
        for i in range(n_datapoints):
            Ks[i] = kernel_rbf(train_x1[i], x_test.numpy(), outputscale, lengthscale)

        # Get variance test points
        Kss = kernel_rbf(x_test.numpy(), x_test.numpy(), outputscale, lengthscale)

        L = np.linalg.cholesky(K)
        v = np.linalg.solve(L, Ks)
        var = Kss - np.matmul(v.T, v)

        model.eval()
        likelihood.eval()
        with gpytorch.settings.fast_pred_var():
            y_preds = likelihood(model(x_test))

        print(f"Predicted variance with gpytorch:{y_preds.variance.item()}")
        print(f"Predicted variance with my kernel:{var}")

4

1 回答 1

0

我发现了错误:

  1. 噪音不是平方的所以它是K += noise * np.eye(n_datapoints)而不是K += noise**2 * np.eye(n_datapoints)
  2. 我忘了在 $$ K** $$ 中添加噪声项,即Kss += noise
于 2021-01-06T10:38:52.467 回答