Repository: cornellius-gp/gpytorch
Pull request #152 (2018-05-31 02:19:05)

Implement white noise kernel
Here's my take on implementing the white noise kernel (#113), and two simple unit tests for it that show it works in at least the simple GP regression case and the 1D KISS-GP regression case.

I changed things a bit so that kernels can return `None`, and `LazyVariables` treat adding `None` as a no-op. If we think this is unsafe or would make debugging more difficult in some cases, we could equivalently make a `ZeroLazyVariable` that behaves the same way.

With that in mind, the white noise kernel returns a `DiagLazyVariable` containing `self.variances` if it believes it is dealing with training data, and returns `None` otherwise. 

There are some changes to `AddedDiagLazyVariable` I'd like to make to make things even more efficient when dealing with this kind of kernel, but that can be a separate PR.

@Balandat 


Patch:
diff --git a/gpytorch/kernels/__init__.py b/gpytorch/kernels/__init__.py
index 83d822d9..24d8cc71 100644
--- a/gpytorch/kernels/__init__.py
+++ b/gpytorch/kernels/__init__.py
@@ -11,7 +11,7 @@ from .linear_kernel import LinearKernel
 from .multiplicative_grid_interpolation_kernel import (
     MultiplicativeGridInterpolationKernel
 )
-
+from .white_noise_kernel import WhiteNoiseKernel
 
 __all__ = [
     Kernel,
@@ -25,4 +25,5 @@ __all__ = [
     InducingPointKernel,
     AdditiveGridInterpolationKernel,
     MultiplicativeGridInterpolationKernel,
+    WhiteNoiseKernel,
 ]
diff --git a/gpytorch/kernels/kernel.py b/gpytorch/kernels/kernel.py
index 77640be2..a579b92d 100644
--- a/gpytorch/kernels/kernel.py
+++ b/gpytorch/kernels/kernel.py
@@ -75,7 +75,7 @@ class AdditiveKernel(Kernel):
         self.kernel_2 = kernel_2
 
     def forward(self, x1, x2):
-        return self.kernel_1(x1, x2) + self.kernel_2(x1, x2)
+        return self.kernel_1(x1, x2).evaluate_kernel() + self.kernel_2(x1, x2).evaluate_kernel()
 
 
 class ProductKernel(Kernel):
@@ -86,4 +86,4 @@ class ProductKernel(Kernel):
         self.kernel_2 = kernel_2
 
     def forward(self, x1, x2):
-        return self.kernel_1(x1, x2) * self.kernel_2(x1, x2)
+        return self.kernel_1(x1, x2).evaluate_kernel() * self.kernel_2(x1, x2).evaluate_kernel()
diff --git a/gpytorch/kernels/white_noise_kernel.py b/gpytorch/kernels/white_noise_kernel.py
new file mode 100644
index 00000000..113d1c46
--- /dev/null
+++ b/gpytorch/kernels/white_noise_kernel.py
@@ -0,0 +1,17 @@
+import torch
+from . import Kernel
+from gpytorch.lazy import DiagLazyVariable, ZeroLazyVariable
+
+
+class WhiteNoiseKernel(Kernel):
+    def __init__(self, variances):
+        super(WhiteNoiseKernel, self).__init__()
+        self.register_buffer("variances", variances)
+
+    def forward(self, x1, x2):
+        if self.training:
+            return DiagLazyVariable(self.variances.unsqueeze(0))
+        elif x1.size(-2) == x2.size(-2) and x1.size(-2) == self.variances.size(-1) and torch.equal(x1, x2):
+            return DiagLazyVariable(self.variances.unsqueeze(0))
+        else:
+            return ZeroLazyVariable(x1.size(-3), x1.size(-2), x1.size(-1), device=x1.device)
diff --git a/gpytorch/lazy/__init__.py b/gpytorch/lazy/__init__.py
index af40b832..8ef785e9 100644
--- a/gpytorch/lazy/__init__.py
+++ b/gpytorch/lazy/__init__.py
@@ -15,6 +15,7 @@ from .root_lazy_variable import RootLazyVariable
 from .sum_lazy_variable import SumLazyVariable
 from .sum_batch_lazy_variable import SumBatchLazyVariable
 from .toeplitz_lazy_variable import ToeplitzLazyVariable
+from .zero_lazy_variable import ZeroLazyVariable
 
 
 __all__ = [
@@ -35,4 +36,5 @@ __all__ = [
     SumLazyVariable,
     SumBatchLazyVariable,
     ToeplitzLazyVariable,
+    ZeroLazyVariable,
 ]
diff --git a/gpytorch/lazy/lazy_evaluated_kernel_variable.py b/gpytorch/lazy/lazy_evaluated_kernel_variable.py
index 7c342a13..a4046129 100644
--- a/gpytorch/lazy/lazy_evaluated_kernel_variable.py
+++ b/gpytorch/lazy/lazy_evaluated_kernel_variable.py
@@ -56,6 +56,7 @@ class LazyEvaluatedKernelVariable(LazyVariable):
             res = super(Kernel, self.kernel).__call__(
                 x1.transpose(-2, -3), x2.transpose(-2, -3)
             )
+
             if isinstance(res, LazyVariable):
                 res = res.evaluate()
             self._cached_kernel_diag = res.transpose(-3, -2).squeeze()
diff --git a/gpytorch/lazy/lazy_variable.py b/gpytorch/lazy/lazy_variable.py
index edb5d38d..8c3a181b 100644
--- a/gpytorch/lazy/lazy_variable.py
+++ b/gpytorch/lazy/lazy_variable.py
@@ -818,13 +818,24 @@ class LazyVariable(object):
 
     def __add__(self, other):
         from .sum_lazy_variable import SumLazyVariable
+        from .zero_lazy_variable import ZeroLazyVariable
+        if isinstance(other, ZeroLazyVariable):
+            return self
 
         return SumLazyVariable(self, other)
 
     def __div__(self, other):
+        from .zero_lazy_variable import ZeroLazyVariable
+        if isinstance(other, ZeroLazyVariable):
+            raise RuntimeError('Attempted to divide by a ZeroLazyVariable (divison by zero)')
+
         return self.mul(1. / other)
 
     def __mul__(self, other):
+        from .zero_lazy_variable import ZeroLazyVariable
+        if isinstance(other, ZeroLazyVariable):
+            return other
+
         return self.mul(other)
 
     def __getitem__(self, index):
diff --git a/gpytorch/lazy/sum_lazy_variable.py b/gpytorch/lazy/sum_lazy_variable.py
index abf07966..f059ff4a 100644
--- a/gpytorch/lazy/sum_lazy_variable.py
+++ b/gpytorch/lazy/sum_lazy_variable.py
@@ -5,6 +5,7 @@ from __future__ import unicode_literals
 
 from .lazy_variable import LazyVariable
 from .non_lazy_variable import NonLazyVariable
+from .zero_lazy_variable import ZeroLazyVariable
 from torch.autograd import Variable
 
 
@@ -89,6 +90,8 @@ class SumLazyVariable(LazyVariable):
         )
 
     def __add__(self, other):
+        if isinstance(other, ZeroLazyVariable):
+            return self
         if isinstance(other, SumLazyVariable):
             return SumLazyVariable(*(list(self.lazy_vars) + list(other.lazy_vars)))
         elif isinstance(other, LazyVariable):
diff --git a/gpytorch/lazy/zero_lazy_variable.py b/gpytorch/lazy/zero_lazy_variable.py
new file mode 100644
index 00000000..25a109e8
--- /dev/null
+++ b/gpytorch/lazy/zero_lazy_variable.py
@@ -0,0 +1,125 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import torch
+from . import LazyVariable
+
+
+class ZeroLazyVariable(LazyVariable):
+    """
+    Special LazyVariable representing zero.
+    """
+    def __init__(self, *sizes, device):
+        super(ZeroLazyVariable, self).__init__(*sizes)
+        self.sizes = list(sizes)
+        self.device = device
+
+    def _matmul(self, rhs):
+        rhs_size_ind = -2 if rhs.ndimension() > 1 else -1
+        if self.size(-1) != rhs.size(rhs_size_ind):
+            raise RuntimeError('Size mismatch, self: {}, rhs: {}'.format(self.size(), rhs.size()))
+        return rhs * 0
+
+    def _t_matmul(self, rhs):
+        rhs_size_ind = -2 if rhs.ndimension() > 1 else -1
+        if self.size(-1) != rhs.size(rhs_size_ind):
+            raise RuntimeError('Size mismatch, self: {}, rhs: {}'.format(self.size(), rhs.size()))
+        return rhs * 0
+
+    def _quad_form_derivative(self, left_vecs, right_vecs):
+        raise RuntimeError('Backwards through a ZeroLazyVariable is not possible')
+
+    def _size(self):
+        return torch.Size(self.sizes)
+
+    def add_diag(self, diag):
+        from .diag_lazy_variable import DiagLazyVariable
+        if diag.size(-1) != self.size(-1):
+            raise RuntimeError('Size mismatch, self: {}, diag {}'.format(self.size(), diag.size()))
+
+        if self.size(-1) != self.size(-2):
+            raise RuntimeError("add_diag only defined for square matrices")
+
+        if self.ndimension() == 3:
+            return DiagLazyVariable(diag.unsqueeze(0).expand(self.size(0), self.size(1)))
+        else:
+            return DiagLazyVariable(diag.expand(self.size(0)))
+
+    def diag(self):
+        size = self.size()
+        if size[-1] != size[-2]:
+            raise RuntimeError("Diag works on square matrices (or batches)")
+
+        if self.ndimension() == 3:
+            return torch.zeros(size[-3], size[-2])
+        else:
+            return torch.zeros(size[-2])
+
+    def evaluate(self):
+        return torch.zeros(*self.sizes, device=self.device)
+
+    def inv_matmul(self, tensor):
+        raise RuntimeError('ZeroLazyVariables are not invertible!')
+
+    def inv_quad(self, tensor):
+        raise RuntimeError('ZeroLazyVariables are not invertible!')
+
+    def inv_quad_log_det(self, inv_quad_rhs=None, log_det=False):
+        raise RuntimeError('ZeroLazyVariables are not invertible!')
+
+    def log_det(self):
+        return torch.log(torch.tensor(0.0))
+
+    def matmul(self, tensor):
+        tensor_size_ind = -2 if tensor.ndimension() > 1 else -1
+        if self.size(-1) != tensor.size(tensor_size_ind):
+            raise RuntimeError('Size mismatch, self: {}, tensor: {}'.format(self.size(), tensor.size()))
+        return tensor * 0
+
+    def mul(self, other):
+        return self
+
+    def mul_batch(self, mul_batch_size=None):
+        return ZeroLazyVariable(*self.sizes[-2:], self.device)
+
+    def root_decomposition(self):
+        raise RuntimeError('ZeroLazyVariables are not positive definite!')
+
+    def root_inv_decomposition(self, initial_vectors=None, test_vectors=None):
+        raise RuntimeError('ZeroLazyVariables are not positive definite!')
+
+    def root_decomposition_size(self):
+        raise RuntimeError('ZeroLazyVariables are not positive definite!')
+
+    def size(self, val=None):
+        """
+        Returns the size of the resulting Variable that the lazy variable represents
+        """
+        size = self._size()
+        if val is not None:
+            return size[val]
+        return size
+
+    def sum_batch(self, sum_batch_size=None):
+        from .sum_batch_lazy_variable import SumBatchLazyVariable
+
+        return SumBatchLazyVariable(self, sum_batch_size=sum_batch_size)
+
+    def transpose(self, dim1, dim2):
+        sizes = self.sizes.copy()
+        tmp = sizes[dim1]
+        sizes[dim1] = sizes[dim2]
+        sizes[dim2] = tmp
+
+        return ZeroLazyVariable(*sizes, device=self.device)
+
+    def __add__(self, other):
+        return other
+
+    def __div__(self, other):
+        return self
+
+    def __mul__(self, other):
+        return self
diff --git a/gpytorch/module.py b/gpytorch/module.py
index 17e3c3bc..71bdbdcb 100644
--- a/gpytorch/module.py
+++ b/gpytorch/module.py
@@ -188,7 +188,6 @@ class Module(nn.Module):
             or isinstance(outputs, LazyVariable)
         ):
             return outputs
-
         for output in outputs:
             if not (
                 isinstance(output, RandomVariable)


Test patch:
diff --git a/test/examples/test_kissgp_white_noise_regression.py b/test/examples/test_kissgp_white_noise_regression.py
new file mode 100644
index 00000000..7c8aae53
--- /dev/null
+++ b/test/examples/test_kissgp_white_noise_regression.py
@@ -0,0 +1,189 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import os
+import math
+import torch
+import unittest
+import gpytorch
+from torch import optim
+from torch.autograd import Variable
+from gpytorch.kernels import RBFKernel, GridInterpolationKernel, WhiteNoiseKernel
+from gpytorch.means import ConstantMean
+from gpytorch.likelihoods import GaussianLikelihood
+from gpytorch.random_variables import GaussianRandomVariable
+
+
+# Simple training data: let's try to learn a sine function,
+# but with KISS-GP let's use 100 training examples.
+def make_data(cuda=False):
+    train_x = Variable(torch.linspace(0, 1, 100))
+    train_y = Variable(torch.sin(train_x.data * (2 * math.pi)))
+    test_x = Variable(torch.linspace(0, 1, 51))
+    test_y = Variable(torch.sin(test_x.data * (2 * math.pi)))
+    if cuda:
+        train_x = train_x.cuda()
+        train_y = train_y.cuda()
+        test_x = test_x.cuda()
+        test_y = test_y.cuda()
+    return train_x, train_y, test_x, test_y
+
+
+class GPRegressionModel(gpytorch.models.ExactGP):
+
+    def __init__(self, train_x, train_y, likelihood):
+        super(GPRegressionModel, self).__init__(train_x, train_y, likelihood)
+        self.mean_module = ConstantMean(constant_bounds=[-1e-5, 1e-5])
+        self.base_covar_module = RBFKernel(log_lengthscale_bounds=(-5, 6))
+        self.grid_covar_module = GridInterpolationKernel(
+            self.base_covar_module, grid_size=50, grid_bounds=[(0, 1)]
+        )
+        self.noise_covar_module = WhiteNoiseKernel(variances=torch.ones(100) * 0.001)
+        self.covar_module = self.grid_covar_module + self.noise_covar_module
+
+    def forward(self, x):
+        mean_x = self.mean_module(x)
+        covar_x = self.covar_module(x)
+        return GaussianRandomVariable(mean_x, covar_x)
+
+
+class TestKISSGPRegression(unittest.TestCase):
+
+    def setUp(self):
+        if (
+            os.getenv("UNLOCK_SEED") is None
+            or os.getenv("UNLOCK_SEED").lower() == "false"
+        ):
+            self.rng_state = torch.get_rng_state()
+            torch.manual_seed(0)
+
+    def tearDown(self):
+        if hasattr(self, "rng_state"):
+            torch.set_rng_state(self.rng_state)
+
+    def test_kissgp_gp_mean_abs_error(self):
+        train_x, train_y, test_x, test_y = make_data()
+        likelihood = GaussianLikelihood()
+        gp_model = GPRegressionModel(train_x.data, train_y.data, likelihood)
+        mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, gp_model)
+
+        # Optimize the model
+        gp_model.train()
+        likelihood.train()
+
+        optimizer = optim.Adam(
+            list(gp_model.parameters()) + list(likelihood.parameters()), lr=0.1
+        )
+        optimizer.n_iter = 0
+        for _ in range(25):
+            optimizer.zero_grad()
+            output = gp_model(train_x)
+            loss = -mll(output, train_y)
+            loss.backward()
+            optimizer.n_iter += 1
+            optimizer.step()
+
+        for param in gp_model.parameters():
+            self.assertTrue(param.grad is not None)
+            self.assertGreater(param.grad.norm().item(), 0)
+        for param in likelihood.parameters():
+            self.assertTrue(param.grad is not None)
+            self.assertGreater(param.grad.norm().item(), 0)
+
+        # Test the model
+        gp_model.eval()
+        likelihood.eval()
+
+        test_preds = likelihood(gp_model(test_x)).mean()
+        mean_abs_error = torch.mean(torch.abs(test_y - test_preds))
+
+        self.assertLess(mean_abs_error.data.squeeze().item(), 0.05)
+
+    def test_kissgp_gp_fast_pred_var(self):
+        with gpytorch.fast_pred_var():
+            train_x, train_y, test_x, test_y = make_data()
+            likelihood = GaussianLikelihood()
+            gp_model = GPRegressionModel(train_x.data, train_y.data, likelihood)
+            mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, gp_model)
+
+            # Optimize the model
+            gp_model.train()
+            likelihood.train()
+
+            optimizer = optim.Adam(
+                list(gp_model.parameters()) + list(likelihood.parameters()), lr=0.1
+            )
+            optimizer.n_iter = 0
+            for _ in range(25):
+                optimizer.zero_grad()
+                output = gp_model(train_x)
+                loss = -mll(output, train_y)
+                loss.backward()
+                optimizer.n_iter += 1
+                optimizer.step()
+
+            for param in gp_model.parameters():
+                self.assertTrue(param.grad is not None)
+                self.assertGreater(param.grad.norm().item(), 0)
+            for param in likelihood.parameters():
+                self.assertTrue(param.grad is not None)
+                self.assertGreater(param.grad.norm().item(), 0)
+
+            # Test the model
+            gp_model.eval()
+            likelihood.eval()
+            # Set the cache
+            test_function_predictions = likelihood(gp_model(train_x))
+
+            # Now bump up the likelihood to something huge
+            # This will make it easy to calculate the variance
+            likelihood.log_noise.data.fill_(3)
+            test_function_predictions = likelihood(gp_model(train_x))
+
+            noise = likelihood.log_noise.exp()
+            var_diff = (test_function_predictions.var() - noise).abs()
+            self.assertLess(torch.max(var_diff.data / noise.data), 0.05)
+
+    def test_kissgp_gp_mean_abs_error_cuda(self):
+        if torch.cuda.is_available():
+            train_x, train_y, test_x, test_y = make_data(cuda=True)
+            likelihood = GaussianLikelihood().cuda()
+            gp_model = GPRegressionModel(train_x.data, train_y.data, likelihood).cuda()
+            mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, gp_model)
+
+            # Optimize the model
+            gp_model.train()
+            likelihood.train()
+
+            optimizer = optim.Adam(
+                list(gp_model.parameters()) + list(likelihood.parameters()), lr=0.1
+            )
+            optimizer.n_iter = 0
+            for _ in range(25):
+                optimizer.zero_grad()
+                output = gp_model(train_x)
+                loss = -mll(output, train_y)
+                loss.backward()
+                optimizer.n_iter += 1
+                optimizer.step()
+
+            for param in gp_model.parameters():
+                self.assertTrue(param.grad is not None)
+                self.assertGreater(param.grad.norm().item(), 0)
+            for param in likelihood.parameters():
+                self.assertTrue(param.grad is not None)
+                self.assertGreater(param.grad.norm().item(), 0)
+
+            # Test the model
+            gp_model.eval()
+            likelihood.eval()
+            test_preds = likelihood(gp_model(test_x)).mean()
+            mean_abs_error = torch.mean(torch.abs(test_y - test_preds))
+
+            self.assertLess(mean_abs_error.data.squeeze().item(), 0.02)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/examples/test_white_noise_regression.py b/test/examples/test_white_noise_regression.py
new file mode 100644
index 00000000..8df09cb1
--- /dev/null
+++ b/test/examples/test_white_noise_regression.py
@@ -0,0 +1,212 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import math
+import torch
+import unittest
+import gpytorch
+from torch import optim
+from torch.autograd import Variable
+from gpytorch.kernels import RBFKernel, WhiteNoiseKernel
+from gpytorch.means import ConstantMean
+from gpytorch.likelihoods import GaussianLikelihood
+from gpytorch.random_variables import GaussianRandomVariable
+
+
+# Simple training data: let's try to learn a sine function
+train_x = Variable(torch.linspace(0, 1, 11))
+train_y = Variable(torch.sin(train_x.data * (2 * math.pi)))
+
+test_x = Variable(torch.linspace(0, 1, 51))
+test_y = Variable(torch.sin(test_x.data * (2 * math.pi)))
+
+
+class ExactGPModel(gpytorch.models.ExactGP):
+
+    def __init__(self, train_inputs, train_targets, likelihood):
+        super(ExactGPModel, self).__init__(train_inputs, train_targets, likelihood)
+        self.mean_module = ConstantMean(constant_bounds=(-1, 1))
+        self.rbf_covar_module = RBFKernel(log_lengthscale_bounds=(-3, 3))
+        self.noise_covar_module = WhiteNoiseKernel(variances=torch.ones(11) * 0.001)
+        self.covar_module = self.rbf_covar_module + self.noise_covar_module
+
+    def forward(self, x):
+        mean_x = self.mean_module(x)
+        covar_x = self.covar_module(x)
+        return GaussianRandomVariable(mean_x, covar_x)
+
+
+class TestSimpleGPRegression(unittest.TestCase):
+
+    def test_posterior_latent_gp_and_likelihood_without_optimization(self):
+        # We're manually going to set the hyperparameters to be ridiculous
+        likelihood = GaussianLikelihood(log_noise_bounds=(-3, 3))
+        gp_model = ExactGPModel(train_x.data, train_y.data, likelihood)
+        # Update bounds to accommodate extreme parameters
+        gp_model.rbf_covar_module.set_bounds(log_lengthscale=(-10, 10))
+        likelihood.set_bounds(log_noise=(-10, 10))
+        # Update parameters
+        gp_model.rbf_covar_module.initialize(log_lengthscale=-10)
+        gp_model.mean_module.initialize(constant=0)
+        likelihood.initialize(log_noise=-10)
+
+        # Compute posterior distribution
+        gp_model.eval()
+        likelihood.eval()
+
+        # Let's see how our model does, conditioned with weird hyperparams
+        # The posterior should fit all the data
+        function_predictions = likelihood(gp_model(train_x))
+
+        self.assertLess(
+            torch.norm(function_predictions.mean().data - train_y.data), 1e-3
+        )
+        self.assertLess(torch.norm(function_predictions.var().data), 5e-3)
+
+        # It shouldn't fit much else though
+        test_function_predictions = gp_model(Variable(torch.Tensor([1.1])))
+
+        self.assertLess(torch.norm(test_function_predictions.mean().data - 0), 1e-4)
+        self.assertLess(torch.norm(test_function_predictions.var().data - 1), 1e-4)
+
+    def test_posterior_latent_gp_and_likelihood_with_optimization(self):
+        # We're manually going to set the hyperparameters to something they shouldn't be
+        likelihood = GaussianLikelihood(log_noise_bounds=(-3, 3))
+        gp_model = ExactGPModel(train_x.data, train_y.data, likelihood)
+        mll = gpytorch.ExactMarginalLogLikelihood(likelihood, gp_model)
+        gp_model.rbf_covar_module.initialize(log_lengthscale=1)
+        gp_model.mean_module.initialize(constant=0)
+        likelihood.initialize(log_noise=1)
+
+        # Find optimal model hyperparameters
+        gp_model.train()
+        likelihood.train()
+        optimizer = optim.Adam(
+            list(gp_model.parameters()) + list(likelihood.parameters()), lr=0.1
+        )
+        optimizer.n_iter = 0
+        for _ in range(50):
+            optimizer.zero_grad()
+            output = gp_model(train_x)
+            loss = -mll(output, train_y)
+            loss.backward()
+            optimizer.n_iter += 1
+            optimizer.step()
+
+        for param in gp_model.parameters():
+            self.assertTrue(param.grad is not None)
+            self.assertGreater(param.grad.norm().item(), 0)
+        for param in likelihood.parameters():
+            self.assertTrue(param.grad is not None)
+            self.assertGreater(param.grad.norm().item(), 0)
+        optimizer.step()
+
+        # Test the model
+        gp_model.eval()
+        likelihood.eval()
+        test_function_predictions = likelihood(gp_model(test_x))
+        mean_abs_error = torch.mean(
+            torch.abs(test_y - test_function_predictions.mean())
+        )
+
+        self.assertLess(mean_abs_error.data.squeeze().item(), 0.05)
+
+    def test_posterior_latent_gp_and_likelihood_fast_pred_var(self):
+        with gpytorch.fast_pred_var():
+            # We're manually going to set the hyperparameters to
+            # something they shouldn't be
+            likelihood = GaussianLikelihood(log_noise_bounds=(-3, 3))
+            gp_model = ExactGPModel(train_x.data, train_y.data, likelihood)
+            mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, gp_model)
+            gp_model.rbf_covar_module.initialize(log_lengthscale=1)
+            gp_model.mean_module.initialize(constant=0)
+            likelihood.initialize(log_noise=1)
+
+            # Find optimal model hyperparameters
+            gp_model.train()
+            likelihood.train()
+            optimizer = optim.Adam(
+                list(gp_model.parameters()) + list(likelihood.parameters()), lr=0.1
+            )
+            optimizer.n_iter = 0
+            for _ in range(50):
+                optimizer.zero_grad()
+                output = gp_model(train_x)
+                loss = -mll(output, train_y)
+                loss.backward()
+                optimizer.n_iter += 1
+                optimizer.step()
+
+            for param in gp_model.parameters():
+                self.assertTrue(param.grad is not None)
+                self.assertGreater(param.grad.norm().item(), 0)
+            for param in likelihood.parameters():
+                self.assertTrue(param.grad is not None)
+                self.assertGreater(param.grad.norm().item(), 0)
+            optimizer.step()
+
+            # Test the model
+            gp_model.eval()
+            likelihood.eval()
+            # Set the cache
+            test_function_predictions = likelihood(gp_model(train_x))
+
+            # Now bump up the likelihood to something huge
+            # This will make it easy to calculate the variance
+            likelihood.log_noise.data.fill_(3)
+            test_function_predictions = likelihood(gp_model(train_x))
+
+            noise = likelihood.log_noise.exp()
+            var_diff = (test_function_predictions.var() - noise).abs()
+
+            self.assertLess(torch.max(var_diff.data / noise.data), 0.05)
+
+    def test_posterior_latent_gp_and_likelihood_with_optimization_cuda(self):
+        if torch.cuda.is_available():
+            # We're manually going to set the hyperparameters to
+            # something they shouldn't be
+            likelihood = GaussianLikelihood(log_noise_bounds=(-3, 3)).cuda()
+            gp_model = ExactGPModel(
+                train_x.data.cuda(), train_y.data.cuda(), likelihood
+            ).cuda()
+            mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, gp_model)
+            gp_model.rbf_covar_module.initialize(log_lengthscale=1)
+            gp_model.mean_module.initialize(constant=0)
+            likelihood.initialize(log_noise=1)
+
+            # Find optimal model hyperparameters
+            gp_model.train()
+            likelihood.train()
+            optimizer = optim.Adam(gp_model.parameters(), lr=0.1)
+            optimizer.n_iter = 0
+            for _ in range(50):
+                optimizer.zero_grad()
+                output = gp_model(train_x.cuda())
+                loss = -mll(output, train_y.cuda())
+                loss.backward()
+                optimizer.n_iter += 1
+                optimizer.step()
+
+            for param in gp_model.parameters():
+                self.assertTrue(param.grad is not None)
+                self.assertGreater(param.grad.norm().item(), 0)
+            for param in likelihood.parameters():
+                self.assertTrue(param.grad is not None)
+                self.assertGreater(param.grad.norm().item(), 0)
+            optimizer.step()
+
+            # Test the model
+            gp_model.eval()
+            likelihood.eval()
+            test_function_predictions = likelihood(gp_model(test_x.cuda()))
+            mean_abs_error = torch.mean(
+                torch.abs(test_y.cuda() - test_function_predictions.mean())
+            )
+
+            self.assertLess(mean_abs_error.data.squeeze().item(), 0.05)
+
+
+if __name__ == "__main__":
+    unittest.main()