Repository: cornellius-gp/gpytorch
Pull request #266 (2018-09-12 21:56:52)

Updates to Kernels
- Introduce `_create_input_grid` method for kernels, which automatically
  creates the grid for pairwise computations
- Update docs for all base kernels
- Backfill unit tests for base kernels
- Refactor SM kernel
- Introduce Cosine kernel

This is a start to fixing #249. It doesn't solve the error, but the next PR/refactor will


Patch:
diff --git a/.travis.yml b/.travis.yml
index 0741be5f..1b132871 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -12,14 +12,14 @@ env:
 
 install:
   - if [[ $TRAVIS_PYTHON_VERSION == 2.7 ]]; then
-            pip install http://download.pytorch.org/whl/cpu/torch-0.4.0-cp27-cp27mu-linux_x86_64.whl;
+            pip install http://download.pytorch.org/whl/cpu/torch-0.4.1-cp27-cp27mu-linux_x86_64.whl;
     else
-            pip install http://download.pytorch.org/whl/cpu/torch-0.4.0-cp36-cp36m-linux_x86_64.whl;
+            pip install http://download.pytorch.org/whl/cpu/torch-0.4.1-cp36-cp36m-linux_x86_64.whl;
     fi
   - python setup.py build develop
 
 script:
-  - python -m unittest discover test/
+  - python -m unittest discover
 
 matrix:
   include:
diff --git a/docs/source/kernels.rst b/docs/source/kernels.rst
index ef973046..aa1f9623 100644
--- a/docs/source/kernels.rst
+++ b/docs/source/kernels.rst
@@ -21,6 +21,12 @@ Kernel
 Standard Kernels
 -----------------------------
 
+:hidden:`CosineKernel`
+~~~~~~~~~~~~~~~~~~~~~~
+
+.. autoclass:: CosineKernel
+   :members:
+
 :hidden:`LinearKernel`
 ~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/gpytorch/kernels/__init__.py b/gpytorch/kernels/__init__.py
index 68a7d6d3..3a0600a3 100644
--- a/gpytorch/kernels/__init__.py
+++ b/gpytorch/kernels/__init__.py
@@ -5,6 +5,7 @@ from __future__ import unicode_literals
 
 from .kernel import Kernel, AdditiveKernel, ProductKernel
 from .additive_grid_interpolation_kernel import AdditiveGridInterpolationKernel
+from .cosine_kernel import CosineKernel
 from .grid_kernel import GridKernel
 from .grid_interpolation_kernel import GridInterpolationKernel
 from .index_kernel import IndexKernel
@@ -23,6 +24,7 @@ __all__ = [
     "Kernel",
     "AdditiveKernel",
     "AdditiveGridInterpolationKernel",
+    "CosineKernel",
     "GridKernel",
     "GridInterpolationKernel",
     "LinearKernel",
diff --git a/gpytorch/kernels/cosine_kernel.py b/gpytorch/kernels/cosine_kernel.py
new file mode 100644
index 00000000..828878a6
--- /dev/null
+++ b/gpytorch/kernels/cosine_kernel.py
@@ -0,0 +1,76 @@
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import math
+import torch
+from .kernel import Kernel
+
+
+class CosineKernel(Kernel):
+    r"""
+    Computes a covariance matrix based on the cosine kernel
+    between inputs :math:`\mathbf{x_1}` and :math:`\mathbf{x_2}`:
+
+    .. math::
+
+       \begin{equation*}
+          k_{\text{Cosine}}(\mathbf{x_1}, \mathbf{x_2}) = \cos \left(
+            \pi \Vert \mathbf{x_1} - \mathbf{x_2} \Vert_2 / p \right)
+       \end{equation*}
+
+    where :math:`p` is the periord length parameter.
+
+    Args:
+        :attr:`batch_size` (int, optional):
+            Set this if you want a separate lengthscale for each
+            batch of input data. It should be `b` if :attr:`x1` is a `b x n x d` tensor. Default: `1`
+        :attr:`active_dims` (tuple of ints, optional):
+            Set this if you want to compute the covariance of only a few input dimensions. The ints
+            corresponds to the indices of the dimensions. Default: `None`.
+        :attr:`log_period_length_prior` (Prior, optional):
+            Set this if you want to apply a prior to the period length parameter.  Default: `None`
+        :attr:`eps` (float):
+            The minimum value that the lengthscale/period length can take
+            (prevents divide by zero errors). Default: `1e-6`.
+
+    Attributes:
+        :attr:`period_length` (Tensor):
+            The period length parameter. Size/shape of parameter depends on the :attr:`batch_size` arguments.
+
+    Example:
+        >>> x = torch.randn(10, 5)
+        >>> # Non-batch: Simple option
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.CosineKernel())
+        >>>
+        >>> batch_x = torch.randn(2, 10, 5)
+        >>> # Batch: Simple option
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.CosineKernel())
+        >>> # Batch: different lengthscale for each batch
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.CosineKernel(batch_size=2))
+        >>> covar = covar_module(x)  # Output: LazyVariable of size (2 x 10 x 10)
+    """
+
+    def __init__(self, active_dims=None, batch_size=1, log_period_length_prior=None, eps=1e-6):
+        super(CosineKernel, self).__init__(has_lengthscale=False, active_dims=active_dims)
+        self.eps = eps
+        self.register_parameter(
+            name="log_period_length",
+            parameter=torch.nn.Parameter(torch.zeros(batch_size, 1, 1)),
+            prior=log_period_length_prior,
+        )
+
+    @property
+    def period_length(self):
+        return self.log_period_length.exp().clamp(self.eps, 1e5)
+
+    def preprocess_data(self, x1, x2):
+        # Override the default behavior:
+        # Here the lengthscale scales the sine, not the inputs directly
+        return x1.div(self.period_length), x2.div(self.period_length)
+
+    def forward(self, x1, x2):
+        diff = torch.norm((x1.unsqueeze(2) - x2.unsqueeze(1)), 2, dim=-1)
+        res = torch.cos(diff.mul(math.pi))
+        return res
diff --git a/gpytorch/kernels/kernel.py b/gpytorch/kernels/kernel.py
index e89874b9..f371ff6a 100644
--- a/gpytorch/kernels/kernel.py
+++ b/gpytorch/kernels/kernel.py
@@ -31,6 +31,52 @@ class Kernel(Module):
         :obj:`torch.tensor` representing the covariance matrix, you may need to call the
         :func:`gpytorch.lazy.LazyTensor.evaluate` method on the output.
 
+    This base :class:`Kernel` class includes a lengthscale parameter
+    :math:`\Theta`, which is used by many common kernel functions.
+    There are a few options for the lengthscale:
+
+    * Default: No lengthscale (i.e. :math:`\Theta` is the identity matrix).
+
+    * Single lengthscale: One lengthscale can be applied to all input dimensions/batches
+      (i.e. :math:`\Theta` is a constant diagonal matrix).
+      This is controlled by setting `has_lengthscale=True`.
+
+    * ARD: Each input dimension gets its own separate lengthscale
+      (i.e. :math:`\Theta` is a non-constant diagonal matrix).
+      This is controlled by the `ard_num_dims` keyword argument (as well has `has_lengthscale=True`).
+
+    In batch-mode (i.e. when :math:`x_1` and :math:`x_2` are batches of input matrices), each
+    batch of data can have its own lengthscale parameter by setting the `batch_size`
+    keyword argument to the appropriate number of batches.
+
+    .. note::
+
+        The :attr:`lengthscale` parameter is parameterized on a log scale to constrain it to be positive.
+        You can set a prior on this parameter using the :attr:`log_lengthscale_prior` argument, but
+        be sure that the prior has :attr:`log_transform=True` set.
+
+    Base Args:
+        :attr:`has_lengthscale` (bool):
+            Set this if the kernel has a lengthscale. Default: `False`.
+        :attr:`ard_num_dims` (int, optional):
+            Set this if you want a separate lengthscale for each input
+            dimension. It should be `d` if :attr:`x1` is a `n x d` matrix.  Default: `None`
+        :attr:`batch_size` (int, optional):
+            Set this if you want a separate lengthscale for each batch of input
+            data. It should be `b` if :attr:`x1` is a `b x n x d` tensor.  Default: `1`
+        :attr:`active_dims` (tuple of ints, optional):
+            Set this if you want to compute the covariance of only a few input dimensions. The ints
+            corresponds to the indices of the dimensions. Default: `None`.
+        :attr:`log_lengthscale_prior` (Prior, optional):
+            Set this if you want to apply a prior to the lengthscale parameter.  Default: `None`
+        :attr:`eps` (float):
+            The minimum value that the lengthscale can take (prevents divide by zero errors). Default: `1e-6`.
+
+    Base Attributes:
+        :attr:`lengthscale` (Tensor):
+            The lengthscale parameter. Size/shape of parameter depends on the
+            :attr:`ard_num_dims` and :attr:`batch_size` arguments.
+
     Example:
         >>> covar_module = gpytorch.kernels.LinearKernel()
         >>> x1 = torch.randn(50, 3)
@@ -42,30 +88,19 @@ class Kernel(Module):
         self,
         has_lengthscale=False,
         ard_num_dims=None,
-        log_lengthscale_prior=None,
-        active_dims=None,
         batch_size=1,
+        active_dims=None,
         log_lengthscale_bounds=None,
+        log_lengthscale_prior=None,
+        eps=1e-6,
     ):
-        """
-        The base Kernel class handles both lengthscales and ARD.
-
-        Args:
-            has_lengthscale (bool): If True, we will register a :obj:`torch.nn.Parameter` named `log_lengthscale`
-            ard_num_dims (int): If not None, the `log_lengthscale` parameter will have this many entries.
-            log_lengthscale_prior (:obj:`gpytorch.priors.Prior`): Prior over the log lengthscale
-            active_dims (list): List of data dimensions to evaluate this Kernel on.
-            batch_size (int): If training or testing multiple GPs simultaneously, this is how many lengthscales to
-                              register.
-            log_lengthscale_bounds (tuple): Deprecated min and max values for the lengthscales. If supplied, this
-                                            now registers a :obj:`gpytorch.priors.SmoothedBoxPrior`
-        """
         super(Kernel, self).__init__()
         if active_dims is not None and not torch.is_tensor(active_dims):
             active_dims = torch.tensor(active_dims, dtype=torch.long)
         self.active_dims = active_dims
         self.ard_num_dims = ard_num_dims
-        self.has_lengthscale = has_lengthscale
+        self.batch_size = batch_size
+        self.__has_lengthscale = has_lengthscale
         if has_lengthscale:
             lengthscale_num_dims = 1 if ard_num_dims is None else ard_num_dims
             log_lengthscale_prior = _bounds_to_prior(prior=log_lengthscale_prior, bounds=log_lengthscale_bounds)
@@ -75,10 +110,14 @@ class Kernel(Module):
                 prior=log_lengthscale_prior,
             )
 
+    @property
+    def has_lengthscale(self):
+        return self.__has_lengthscale
+
     @property
     def lengthscale(self):
-        if "log_lengthscale" in self.named_parameters().keys():
-            return self.log_lengthscale.exp()
+        if self.has_lengthscale:
+            return self.log_lengthscale.exp().clamp(self.eps, 1e5)
         else:
             return None
 
@@ -103,28 +142,77 @@ class Kernel(Module):
 
     @abstractmethod
     def forward(self, x1, x2, **params):
+        """
+        Computes the covariance between x1 and x2.
+        This method should be imlemented by all Kernel subclasses.
+
+        .. note::
+
+            All non-compositional kernels should use the :meth:`gpytorch.kernels.Kernel._create_input_grid`
+            method to create a meshgrid between x1 and x2 (if necessary).
+
+            Do not manually create the grid - this is inefficient and will cause erroneous behavior in certain
+            evaluation modes.
+
+        Args:
+            :attr:`x1` (Tensor `n x d` or `b x n x d`)
+            :attr:`x2` (Tensor `m x d` or `b x m x d`) - for diag mode, these must be the same inputs
+
+        Returns:
+            :class:`Tensor` or :class:`gpytorch.lazy.LazyTensor`.
+            The exact size depends on the kernel's evaluation mode:
+
+            * `full_covar`: `n x m` or `b x n x m`
+            * `full_covar` with `dim_groups=k`: `k x n x m` or `b x k x n x m`
+            * `diag`: `n` or `b x n`
+            * `diag` with `dim_groups=k`: `k x n` or `b x k x n`
+        """
         raise NotImplementedError()
 
-    def __call__(self, x1_, x2_=None, **params):
-        x1, x2 = x1_, x2_
+    def _create_input_grid(self, x1, x2):
+        """
+        This is a helper method for creating a grid of the kernel's inputs.
+        Use this helper rather than maually creating a meshgrid.
+
+        The grid dimensions depend on the kernel's evaluation mode.
+
+        Args:
+            :attr:`x1` (Tensor `n x d` or `b x n x d`)
+            :attr:`x2` (Tensor `m x d` or `b x m x d`) - for diag mode, these must be the same inputs
+
+        Returns:
+            (:class:`Tensor`, :class:`Tensor) corresponding to the gridded `x1` and `x2`.
+            The shape depends on the kernel's mode
+
+            * `full_covar`: (`b x n x 1 x d` and `b x 1 x m x d`)
+            * `full_covar` with `dim_groups=k`: (`b x k x n x 1 x (d//k)` and `b x k x 1 x m x (d//k)`)
+            * `diag`: (`b x n x d` and `b x n x d`)
+            * `diag` with `dim_groups=k`: (`b x k x n x (d//k)` and `b x k x n x (d//k)`)
+        """
+        return x1.unsqueeze(-2), x2.unsqueeze(-3)
+
+    def __call__(self, x1, x2=None, **params):
+        x1_, x2_ = x1, x2
 
+        # Select the active dimensions
         if self.active_dims is not None:
-            x1 = x1_.index_select(-1, self.active_dims)
+            x1_ = x1_.index_select(-1, self.active_dims)
             if x2_ is not None:
-                x2 = x2_.index_select(-1, self.active_dims)
+                x2_ = x2_.index_select(-1, self.active_dims)
 
-        if x2 is None:
-            x2 = x1
+        # Give x1_ and x2_ a last dimension, if necessary
+        if x1_.ndimension() == 1:
+            x1_ = x1_.unsqueeze(1)
+        if x2_ is not None:
+            if x2_.ndimension() == 1:
+                x2_ = x2_.unsqueeze(1)
+            if not x1_.size(-1) == x2_.size(-1):
+                raise RuntimeError("x1_ and x2_ must have the same number of dimensions!")
 
-        # Give x1 and x2 a last dimension, if necessary
-        if x1.ndimension() == 1:
-            x1 = x1.unsqueeze(1)
-        if x2.ndimension() == 1:
-            x2 = x2.unsqueeze(1)
-        if not x1.size(-1) == x2.size(-1):
-            raise RuntimeError("x1 and x2 must have the same number of dimensions!")
+        if x2_ is None:
+            x2_ = x1_
 
-        return LazyEvaluatedKernelTensor(self, x1, x2)
+        return LazyEvaluatedKernelTensor(self, x1_, x2_, **params)
 
     def __add__(self, other):
         return AdditiveKernel(self, other)
diff --git a/gpytorch/kernels/linear_kernel.py b/gpytorch/kernels/linear_kernel.py
index 684d0c75..794df096 100644
--- a/gpytorch/kernels/linear_kernel.py
+++ b/gpytorch/kernels/linear_kernel.py
@@ -10,14 +10,42 @@ from gpytorch.priors._compatibility import _bounds_to_prior
 
 
 class LinearKernel(Kernel):
-    """
-    An implementation of the linear kernel :math:`k(x, z) =(x-offset)(z-offset)' + variance`.
+    r"""
+    Computes a covariance matrix based on the Linear kernel
+    between inputs :math:`\mathbf{x_1}` and :math:`\mathbf{x_2}`:
+
+    .. math::
+        \begin{equation*}
+            k_\text{Linear}(\mathbf{x_1}, \mathbf{x_2}) = (\mathbf{x_1} - \mathbf{o})^\top
+            (\mathbf{x_2} - \mathbf{o}) + v.
+        \end{equation*}
+
+    where
+
+    * :math:`\mathbf o` is an :attr:`offset` parameter.
+    * :math:`v` is a :attr:`variance` parameter.
+
+
+    .. note::
+
+        To implement this efficiently, we use a :obj:`gpytorch.lazy.RootLazyTensor` during training and a
+        :class:`gpytorch.lazy.MatmulLazyTensor` during test. These lazy tensors represent matrices of the form
+        :math:`K = XX^{\top}` and :math:`K = XZ^{\top}`. This makes inference
+        efficient because a matrix-vector product :math:`Kv` can be computed as
+        :math:`Kv=X(X^{\top}v)`, where the base multiply :math:`Xv` takes only
+        :math:`O(nd)` time and space.
 
-    To implement this efficiently, we use a :obj:`gpytorch.lazy.RootLazyTensor` during training and a
-    :math:`gpytorch.lazy.MatmulLazyTensor` during test. These lazy tensors represent matrices of the form
-    :math:`K = XX^{\top}` and :math:`K = XZ^{\top}`. This makes inference efficient because a matrix-vector product
-    :math:`Kv` can be computed as :math:`Kv=X(X^{\top}v)`, where the base multiply :math:`Xv` takes only :math:`O(nd)`
-    time and space.
+    Args:
+        :attr:`num_dimensions` (int):
+            Number of data dimensions to expect. This
+            is necessary to create the offset parameter.
+        :attr:`variance_prior` (:class:`gpytorch.priors.Prior`):
+            Prior over the variance parameter (default `None`).
+        :attr:`offset_prior` (:class:`gpytorch.priors.Prior`):
+            Prior over the offset parameter (default `None`).
+        :attr:`active_dims` (list):
+            List of data dimensions to operate on.
+            `len(active_dims)` should equal `num_dimensions`.
     """
 
     def __init__(
@@ -29,17 +57,6 @@ class LinearKernel(Kernel):
         variance_bounds=None,
         offset_bounds=None,
     ):
-        """
-        Args:
-            num_dimensions (int): Number of data dimensions to expect. This is necessary to create the offset parameter.
-            variance_prior (:obj:`gpytorch.priors.Prior`): Prior over the variance parameter (default `None`).
-            offset_prior (:obj:`gpytorch.priors.Prior`): Prior over the offset parameter (default `None`).
-            active_dims (list): List of data dimensions to operate on. `len(active_dims)` should equal `num_dimensions`.
-            variance_bounds (tuple, deprecated): Min and max value for the variance parameter. Deprecated, and now
-                                                 creates a :obj:`gpytorch.priors.SmoothedBoxPrior`.
-            offset_bounds (tuple, deprecated): Min and max value for the offset parameter. Deprecated, and now creates a
-                                                :obj:'gpytorch.priors.SmoothedBoxPrior'.
-        """
         super(LinearKernel, self).__init__(active_dims=active_dims)
         variance_prior = _bounds_to_prior(prior=variance_prior, bounds=variance_bounds, log_transform=False)
         self.register_parameter(name="variance", parameter=torch.nn.Parameter(torch.zeros(1)), prior=variance_prior)
diff --git a/gpytorch/kernels/matern_kernel.py b/gpytorch/kernels/matern_kernel.py
index 141989d9..388c1444 100644
--- a/gpytorch/kernels/matern_kernel.py
+++ b/gpytorch/kernels/matern_kernel.py
@@ -9,14 +9,82 @@ from .kernel import Kernel
 
 
 class MaternKernel(Kernel):
+    r"""
+    Computes a covariance matrix based on the Matern kernel
+    between inputs :math:`\mathbf{x_1}` and :math:`\mathbf{x_2}`:
+
+    .. math::
+
+       \begin{equation*}
+          k_{\text{Matern}}(\mathbf{x_1}, \mathbf{x_2}) = \frac{2^{1 - \nu}}{\Gamma(\nu)}
+          \left( \sqrt{2 \nu} d \right) K_\nu \left( \sqrt{2 \nu} d \right)
+       \end{equation*}
+
+    where
+
+    * :math:`d = (\mathbf{x_1} - \mathbf{x_2})^\top \Theta^{-1} (\mathbf{x_1} - \mathbf{x_2})`
+      is the distance between
+      :math:`x_1` and :math:`x_2` scaled by the :attr:`lengthscale` parameter :math:`\Theta`.
+    * :math:`\nu` is a smoothness parameter (takes values 1/2, 3/2, or 5/2). Smaller values are less smooth.
+    * :math:`K_\nu` is a modified Bessel function.
+
+    There are a few options for the lengthscale parameter :math:`\Theta`:
+    See :class:`gpytorch.kernels.Kernel` for descriptions of the lengthscale options.
+
+    .. note::
+
+        This kernel does not have an `outputscale` parameter. To add a scaling parameter,
+        decorate this kernel with a :class:`gpytorch.kernels.ScaleKernel`.
+
+    Args:
+        :attr:`nu` (float):
+            The smoothness parameter: either 1/2, 3/2, or 5/2.
+        :attr:`ard_num_dims` (int, optional):
+            Set this if you want a separate lengthscale for each
+            input dimension. It should be `d` if :attr:`x1` is a `n x d` matrix. Default: `None`
+        :attr:`batch_size` (int, optional):
+            Set this if you want a separate lengthscale for each
+            batch of input data. It should be `b` if :attr:`x1` is a `b x n x d` tensor. Default: `1`
+        :attr:`active_dims` (tuple of ints, optional):
+            Set this if you want to
+            compute the covariance of only a few input dimensions. The ints
+            corresponds to the indices of the dimensions. Default: `None`.
+        :attr:`log_lengthscale_prior` (Prior, optional):
+            Set this if you want
+            to apply a prior to the lengthscale parameter.  Default: `None`
+        :attr:`eps` (float):
+            The minimum value that the lengthscale can take
+            (prevents divide by zero errors). Default: `1e-6`.
+
+    Attributes:
+        :attr:`lengthscale` (Tensor):
+            The lengthscale parameter. Size/shape of parameter depends on the
+            :attr:`ard_num_dims` and :attr:`batch_size` arguments.
+
+    Example:
+        >>> x = torch.randn(10, 5)
+        >>> # Non-batch: Simple option
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.MaternKernel(nu=0.5))
+        >>> # Non-batch: ARD (different lengthscale for each input dimension)
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.MaternKernel(nu=0.5, ard_num_dims=5))
+        >>> covar = covar_module(x)  # Output: LazyVariable of size (10 x 10)
+        >>>
+        >>> batch_x = torch.randn(2, 10, 5)
+        >>> # Batch: Simple option
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.MaternKernel(nu=0.5))
+        >>> # Batch: different lengthscale for each batch
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.MaternKernel(nu=0.5, batch_size=2))
+        >>> covar = covar_module(x)  # Output: LazyVariable of size (2 x 10 x 10)
+    """
+
     def __init__(
         self,
         nu=2.5,
         ard_num_dims=None,
-        log_lengthscale_prior=None,
-        active_dims=None,
-        eps=1e-8,
         batch_size=1,
+        active_dims=None,
+        eps=1e-6,
+        log_lengthscale_prior=None,
         log_lengthscale_bounds=None,
     ):
         if nu not in {0.5, 1.5, 2.5}:
@@ -33,10 +101,9 @@ class MaternKernel(Kernel):
         self.eps = eps
 
     def forward_diag(self, x1, x2):
-        lengthscale = self.log_lengthscale.exp()
-        mean = x1.mean(1).mean(0)
-        x1_normed = (x1 - mean.unsqueeze(0).unsqueeze(1)).div(lengthscale)
-        x2_normed = (x2 - mean.unsqueeze(0).unsqueeze(1)).div(lengthscale)
+        mean = x1.mean(1, keepdim=True).mean(0, keepdim=True)
+        x1_normed = x1 - mean.unsqueeze(0).unsqueeze(1)
+        x2_normed = x2 - mean.unsqueeze(0).unsqueeze(1)
 
         diff = x1_normed - x2_normed
         distance_over_rho = diff.pow_(2).sum(-1).sqrt()
@@ -52,23 +119,18 @@ class MaternKernel(Kernel):
         return constant_component * exp_component
 
     def forward(self, x1, x2):
-        lengthscale = self.log_lengthscale.exp()
-        mean = x1.mean(1).mean(0)
-        x1_normed = (x1 - mean.unsqueeze(0).unsqueeze(1)).div(lengthscale)
-        x2_normed = (x2 - mean.unsqueeze(0).unsqueeze(1)).div(lengthscale)
+        mean = x1.view(-1, 1, *list(x1.size())[2:]).mean(0, keepdim=True)
+        x1_, x2_ = self._create_input_grid(x1 - mean, x2 - mean)
+        x1_ = x1_.div(self.lengthscale)
+        x2_ = x2_.div(self.lengthscale)
 
-        x1_squared = x1_normed.norm(2, -1).pow(2)
-        x2_squared = x2_normed.norm(2, -1).pow(2)
-        x1_t_x_2 = torch.matmul(x1_normed, x2_normed.transpose(-1, -2))
-
-        distance_over_rho = x1_squared.unsqueeze(-1) + x2_squared.unsqueeze(-2) - x1_t_x_2.mul(2)
-        distance_over_rho = distance_over_rho.clamp(self.eps, 1e10).sqrt()
-        exp_component = torch.exp(-math.sqrt(self.nu * 2) * distance_over_rho)
+        distance = (x1_ - x2_).norm(2, dim=-1)
+        exp_component = torch.exp(-math.sqrt(self.nu * 2) * distance)
 
         if self.nu == 0.5:
             constant_component = 1
         elif self.nu == 1.5:
-            constant_component = (math.sqrt(3) * distance_over_rho).add(1)
+            constant_component = (math.sqrt(3) * distance).add(1)
         elif self.nu == 2.5:
-            constant_component = (math.sqrt(5) * distance_over_rho).add(1).add(5. / 3. * distance_over_rho ** 2)
+            constant_component = (math.sqrt(5) * distance).add(1).add(5. / 3. * distance ** 2)
         return constant_component * exp_component
diff --git a/gpytorch/kernels/multitask_kernel.py b/gpytorch/kernels/multitask_kernel.py
index 74a64156..7955435f 100644
--- a/gpytorch/kernels/multitask_kernel.py
+++ b/gpytorch/kernels/multitask_kernel.py
@@ -60,7 +60,7 @@ class MultitaskKernel(Kernel):
 
     def forward(self, x1, x2):
         covar_i = self.task_covar_module.covar_matrix
-        covar_x = self.data_covar_module.forward(x1, x2)
+        covar_x = self.data_covar_module(x1, x2)
         if covar_x.size(0) == 1:
             covar_x = covar_x[0]
         if not isinstance(covar_x, LazyTensor):
diff --git a/gpytorch/kernels/periodic_kernel.py b/gpytorch/kernels/periodic_kernel.py
index fbda08bc..9219a11b 100644
--- a/gpytorch/kernels/periodic_kernel.py
+++ b/gpytorch/kernels/periodic_kernel.py
@@ -5,17 +5,82 @@ from __future__ import unicode_literals
 
 import math
 import torch
-from gpytorch.kernels import Kernel
-from gpytorch.priors._compatibility import _bounds_to_prior
+from .kernel import Kernel
+from ..priors._compatibility import _bounds_to_prior
 
 
 class PeriodicKernel(Kernel):
+    r""" Computes a covariance matrix based on the periodic kernel
+    between inputs :math:`\mathbf{x_1}` and :math:`\mathbf{x_2}`:
+
+    .. math::
+
+       \begin{equation*}
+          k_{\text{Periodic}}(\mathbf{x_1}, \mathbf{x_2}) = \exp \left(
+            \frac{2 \sin^2 \left( \pi \Vert \mathbf{x_1} - \mathbf{x_2} \Vert_1 / p \right) }
+            { \ell^2 } \right)
+       \end{equation*}
+
+    where
+
+    * :math:`p` is the periord length parameter.
+    * :math:`\ell` is a lengthscale parameter.
+
+    .. note::
+
+        This kernel does not have an `outputscale` parameter. To add a scaling parameter,
+        decorate this kernel with a :class:`gpytorch.kernels.ScaleKernel`.
+
+    .. note::
+
+        This kernel does not have an ARD lengthscale option.
+
+    Args:
+        :attr:`batch_size` (int, optional):
+            Set this if you want a separate lengthscale for each
+            batch of input data. It should be `b` if :attr:`x1` is a `b x n x d` tensor. Default: `1`
+        :attr:`active_dims` (tuple of ints, optional):
+            Set this if you want to
+            compute the covariance of only a few input dimensions. The ints
+            corresponds to the indices of the dimensions. Default: `None`.
+        :attr:`log_period_length_prior` (Prior, optional):
+            Set this if you want
+            to apply a prior to the period length parameter.  Default: `None`
+        :attr:`log_lengthscale_prior` (Prior, optional):
+            Set this if you want
+            to apply a prior to the lengthscale parameter.  Default: `None`
+        :attr:`eps` (float):
+            The minimum value that the lengthscale/period length can take
+            (prevents divide by zero errors). Default: `1e-6`.
+
+    Attributes:
+        :attr:`lengthscale` (Tensor):
+            The lengthscale parameter. Size/shape of parameter depends on the
+            :attr:`batch_size` arguments.
+        :attr:`period_length` (Tensor):
+            The period length parameter. Size/shape of parameter depends on the
+            :attr:`batch_size` arguments.
+
+    Example:
+        >>> x = torch.randn(10, 5)
+        >>> # Non-batch: Simple option
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.PeriodicKernel())
+        >>>
+        >>> batch_x = torch.randn(2, 10, 5)
+        >>> # Batch: Simple option
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.PeriodicKernel())
+        >>> # Batch: different lengthscale for each batch
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.PeriodicKernel(batch_size=2))
+        >>> covar = covar_module(x)  # Output: LazyVariable of size (2 x 10 x 10)
+    """
+
     def __init__(
         self,
+        active_dims=None,
+        batch_size=1,
+        eps=1e-6,
         log_lengthscale_prior=None,
         log_period_length_prior=None,
-        eps=1e-5,
-        active_dims=None,
         log_lengthscale_bounds=None,
         log_period_length_bounds=None,
     ):
@@ -28,12 +93,20 @@ class PeriodicKernel(Kernel):
         )
         self.eps = eps
         self.register_parameter(
-            name="log_period_length", parameter=torch.nn.Parameter(torch.zeros(1, 1, 1)), prior=log_period_length_prior
+            name="log_period_length",
+            parameter=torch.nn.Parameter(torch.zeros(batch_size, 1, 1)),
+            prior=log_period_length_prior,
         )
 
+    @property
+    def period_length(self):
+        return self.log_period_length.exp().clamp(self.eps, 1e5)
+
     def forward(self, x1, x2):
-        lengthscale = (self.log_lengthscale.exp() + self.eps).sqrt_()
-        period_length = (self.log_period_length.exp() + self.eps).sqrt_()
-        diff = torch.sum((x1.unsqueeze(2) - x2.unsqueeze(1)).abs(), -1)
-        res = -2 * torch.sin(math.pi * diff / period_length).pow(2) / lengthscale
-        return res.exp()
+        x1_, x2_ = self._create_input_grid(x1, x2)
+        x1_ = x1_.div(self.period_length)
+        x2_ = x2_.div(self.period_length)
+
+        diff = torch.sum((x1_ - x2_).abs(), -1)
+        res = torch.sin(diff.mul(math.pi)).pow(2).mul(-2 / self.lengthscale).exp_()
+        return res
diff --git a/gpytorch/kernels/rbf_kernel.py b/gpytorch/kernels/rbf_kernel.py
index 7a5e1438..bd21e70b 100644
--- a/gpytorch/kernels/rbf_kernel.py
+++ b/gpytorch/kernels/rbf_kernel.py
@@ -10,33 +10,17 @@ from .kernel import Kernel
 class RBFKernel(Kernel):
     r"""
     Computes a covariance matrix based on the RBF (squared exponential) kernel
-    between inputs :math:`x_1` and :math:`x_2`:
+    between inputs :math:`\mathbf{x_1}` and :math:`\mathbf{x_2}`:
 
     .. math::
 
        \begin{equation*}
-          K_{\text{RBF}} = \exp \left( -\frac{1}{2} (x_1 - x_2)^\top \Theta^{-1} (x_1 - x_2) \right)
+          k_{\text{RBF}}(\mathbf{x_1}, \mathbf{x_2}) = \exp \left( -\frac{1}{2}
+          (\mathbf{x_1} - \mathbf{x_2})^\top \Theta^{-1} (\mathbf{x_1} - \mathbf{x_2}) \right)
        \end{equation*}
 
     where :math:`\Theta` is a :attr:`lengthscale` parameter.
-    There are a few options for the lengthscale:
-
-    * Default: One lengthscale can be applied to all input dimensions/batches
-      (i.e. :math:`\Theta` is a constant diagonal matrix).
-
-    * ARD: Each input dimension gets its own separate lengthscale
-      (i.e. :math:`\Theta` is a non-constant diagonal matrix).
-      This is controlled by the `ard_num_dims` keyword argument.
-
-    In batch-mode (i.e. when :math:`x_1` and :math:`x_2` are batches of input matrices), each
-    batch of data can have its own lengthscale parameter by setting the `batch_size`
-    keyword argument to the appropriate number of batches.
-
-    .. note::
-
-        The :attr:`lengthscale` parameter is parameterized on a log scale to constrain it to be positive.
-        You can set a prior on this parameter using the :attr:`log_lengthscale_prior` argument, but
-        be sure that the prior has :attr:`log_transform=True` set.
+    See :class:`gpytorch.kernels.Kernel` for descriptions of the lengthscale options.
 
     .. note::
 
@@ -44,35 +28,41 @@ class RBFKernel(Kernel):
         decorate this kernel with a :class:`gpytorch.kernels.ScaleKernel`.
 
     Args:
-        :attr:`ard_num_dims` (int, optional): Set this if you want a separate lengthscale for each
+        :attr:`ard_num_dims` (int, optional):
+            Set this if you want a separate lengthscale for each
             input dimension. It should be `d` if :attr:`x1` is a `n x d` matrix. Default: `None`
-        :attr:`batch_size` (int, optional): Set this if you want a separate lengthscale for each
+        :attr:`batch_size` (int, optional):
+            Set this if you want a separate lengthscale for each
             batch of input data. It should be `b` if :attr:`x1` is a `b x n x d` tensor. Default: `1`
-        :attr:`active_dims` (tuple of ints, optional): Set this if you want to
+        :attr:`active_dims` (tuple of ints, optional):
+            Set this if you want to
             compute the covariance of only a few input dimensions. The ints
             corresponds to the indices of the dimensions. Default: `None`.
-        :attr:`log_lengthscale_prior` (Prior, optional): Set this if you want
+        :attr:`log_lengthscale_prior` (Prior, optional):
+            Set this if you want
             to apply a prior to the lengthscale parameter.  Default: `None`
-        :attr:`eps` (float): The minimum value that the lengthscale can take
+        :attr:`eps` (float):
+            The minimum value that the lengthscale can take
             (prevents divide by zero errors). Default: `1e-6`.
 
     Attributes:
-        :attr:`lengthscale` (Tensor): The lengthscale parameter. Size/shape of parameter depends on the
+        :attr:`lengthscale` (Tensor):
+            The lengthscale parameter. Size/shape of parameter depends on the
             :attr:`ard_num_dims` and :attr:`batch_size` arguments.
 
     Example:
         >>> x = torch.randn(10, 5)
         >>> # Non-batch: Simple option
-        >>> covar_module = gpytorch.kernels.RBFKernel()
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())
         >>> # Non-batch: ARD (different lengthscale for each input dimension)
-        >>> covar_module = gpytorch.kernels.RBFKernel(ard_num_dims=5)
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel(ard_num_dims=5))
         >>> covar = covar_module(x)  # Output: LazyTensor of size (10 x 10)
         >>>
         >>> batch_x = torch.randn(2, 10, 5)
         >>> # Batch: Simple option
-        >>> covar_module = gpytorch.kernels.RBFKernel()
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())
         >>> # Batch: different lengthscale for each batch
-        >>> covar_module = gpytorch.kernels.RBFKernel(batch_size=2)
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel(batch_size=2))
         >>> covar = covar_module(x)  # Output: LazyTensor of size (2 x 10 x 10)
     """
 
@@ -102,6 +92,9 @@ class RBFKernel(Kernel):
         return res
 
     def forward(self, x1, x2):
-        lengthscales = self.log_lengthscale.exp().mul(math.sqrt(2)).clamp(self.eps, 1e5)
-        diff = (x1.unsqueeze(2) - x2.unsqueeze(1)).div_(lengthscales.unsqueeze(1))
-        return diff.pow_(2).sum(-1).mul_(-1).exp_()
+        x1_, x2_ = self._create_input_grid(x1, x2)
+        x1_ = x1_.div(self.lengthscale)
+        x2_ = x2_.div(self.lengthscale)
+
+        diff = (x1_ - x2_).norm(2, dim=-1)
+        return diff.pow(2).div_(-2).exp_()
diff --git a/gpytorch/kernels/spectral_mixture_kernel.py b/gpytorch/kernels/spectral_mixture_kernel.py
index 68976b18..127a7823 100644
--- a/gpytorch/kernels/spectral_mixture_kernel.py
+++ b/gpytorch/kernels/spectral_mixture_kernel.py
@@ -3,41 +3,133 @@ from __future__ import division
 from __future__ import print_function
 from __future__ import unicode_literals
 
+import warnings
 import logging
 import math
 import torch
-from gpytorch.kernels import Kernel
+from .kernel import Kernel
 
 logger = logging.getLogger()
 
 
 class SpectralMixtureKernel(Kernel):
+    r"""
+    Computes a covariance matrix based on the Spectral Mixture Kernel
+    between inputs :math:`mathbf{x_1}` and :math:`mathbf{x_2}`:
+    It was proposed in `Gaussian Process Kernels for Pattern Discovery and Extrapolation`_.
+
+    .. note::
+
+        Unlike other kernels,
+        * :attr:`ard_num_dums` **must equal** the number of dimensions of the data
+        * :attr:`batch_size` **must equal** the batch size of the data (1 if the data is not batched)
+        * This kernel should not be combined with a :class:`gpytorch.kernels.ScaleKernel`.
+
+    Args:
+        :attr:`num_mixtures` (int, optional):
+            The number of components in the mixture.
+        :attr:`ard_num_dims` (int, optional):
+            Set this to match the dimensionality of the input.
+            It should be `d` if :attr:`x1` is a `n x d` matrix. Default: `1`
+        :attr:`batch_size` (int, optional):
+            Set this if the data is
+            batch of input data. It should be `b` if :attr:`x1` is a `b x n x d` tensor. Default: `1`
+        :attr:`active_dims` (tuple of ints, optional):
+            Set this if you want to
+            compute the covariance of only a few input dimensions. The ints
+            corresponds to the indices of the dimensions. Default: `None`.
+        :attr:`eps` (float):
+            The minimum value that the lengthscale can take
+            (prevents divide by zero errors). Default: `1e-6`.
+
+    Attributes:
+        :attr:`mixture_lengthscale` (Tensor):
+            The lengthscale parameter. Given `k` mixture components,
+            and `b x n x d` data, this will be of size `b x k x 1 x d`.
+        :attr:`mixture_means` (Tensor):
+            The mixture mean parameters (`b x k x 1 x d`).
+        :attr:`mixture_weights` (Tensor):
+            The mixture weight parameters (`b x k`).
+
+    Example:
+        >>> x = torch.randn(10, 5)
+        >>> # Non-batch: Simple option
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())
+        >>> # Non-batch: ARD (different lengthscale for each input dimension)
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel(ard_num_dims=5))
+        >>> covar = covar_module(x)  # Output: LazyVariable of size (10 x 10)
+        >>>
+        >>> batch_x = torch.randn(2, 10, 5)
+        >>> # Batch: Simple option
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())
+        >>> # Batch: different lengthscale for each batch
+        >>> covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel(batch_size=2))
+        >>> covar = covar_module(x)  # Output: LazyVariable of size (2 x 10 x 10)
+
+
+    .. Gaussian Process Kernels for Pattern Discovery and Extrapolation:
+        https://arxiv.org/pdf/1302.4245.pdf
+    """
+    # TODO: add equation to docs
+
     def __init__(
         self,
-        n_mixtures,
-        n_dims=1,
-        log_mixture_weight_prior=None,
-        log_mixture_mean_prior=None,
-        log_mixture_scale_prior=None,
+        num_mixtures=None,
+        ard_num_dims=1,
+        batch_size=1,
         active_dims=None,
+        eps=1e-6,
+        log_mixture_scales_prior=None,
+        log_mixture_means_prior=None,
+        log_mixture_weights_prior=None,
+        n_mixtures=None,
+        n_dims=None,
     ):
-        self.n_mixtures = n_mixtures
-        self.n_dims = n_dims
+        if n_mixtures is not None:
+            warnings.warn("n_mixtures is deprecated. Use num_mixtures instead.", DeprecationWarning)
+            num_mixtures = n_mixtures
+        if num_mixtures is None:
+            raise RuntimeError("num_mixtures is a required argument")
+        if n_dims is not None:
+            warnings.warn("n_dims is deprecated. Use ard_num_dims instead.", DeprecationWarning)
+            ard_num_dims = n_dims
         if (
-            log_mixture_mean_prior is not None
-            or log_mixture_scale_prior is not None
-            or log_mixture_weight_prior is not None
+            log_mixture_means_prior is not None
+            or log_mixture_scales_prior is not None
+            or log_mixture_weights_prior is not None
         ):
             logger.warning("Priors not implemented for SpectralMixtureKernel")
 
+        # This kernel does not use the default lengthscale
         super(SpectralMixtureKernel, self).__init__(active_dims=active_dims)
-        self.register_parameter(name="log_mixture_weights", parameter=torch.nn.Parameter(torch.zeros(self.n_mixtures)))
+        self.num_mixtures = num_mixtures
+        self.batch_size = batch_size
+        self.ard_num_dims = ard_num_dims
+        self.eps = eps
+
         self.register_parameter(
-            name="log_mixture_means", parameter=torch.nn.Parameter(torch.zeros(self.n_mixtures, self.n_dims))
+            name="log_mixture_weights", parameter=torch.nn.Parameter(torch.zeros(self.batch_size, self.num_mixtures))
         )
         self.register_parameter(
-            name="log_mixture_scales", parameter=torch.nn.Parameter(torch.zeros(self.n_mixtures, self.n_dims))
+            name="log_mixture_means",
+            parameter=torch.nn.Parameter(torch.zeros(self.batch_size, self.num_mixtures, 1, self.ard_num_dims)),
         )
+        self.register_parameter(
+            name="log_mixture_scales",
+            parameter=torch.nn.Parameter(torch.zeros(self.batch_size, self.num_mixtures, 1, self.ard_num_dims)),
+        )
+
+    @property
+    def mixture_scales(self):
+        return self.log_mixture_scales.exp().clamp(self.eps, 1e5)
+
+    @property
+    def mixture_means(self):
+        return self.log_mixture_means.exp().clamp(self.eps, 1e5)
+
+    @property
+    def mixture_weights(self):
+        return self.log_mixture_weights.exp().clamp(self.eps, 1e5)
 
     def initialize_from_data(self, train_x, train_y, **kwargs):
         if not torch.is_tensor(train_x) or not torch.is_tensor(train_y):
@@ -50,8 +142,8 @@ class SpectralMixtureKernel(Kernel):
         train_x_sort = train_x.sort(1)[0]
         max_dist = train_x_sort[:, -1, :] - train_x_sort[:, 0, :]
         min_dist_sort = (train_x_sort[:, 1:, :] - train_x_sort[:, :-1, :]).squeeze(0)
-        min_dist = torch.zeros(1, self.n_dims)
-        for ind in range(self.n_dims):
+        min_dist = torch.zeros(1, self.ard_num_dims)
+        for ind in range(self.ard_num_dims):
             min_dist[:, ind] = min_dist_sort[(torch.nonzero(min_dist_sort[:, ind]))[0], ind]
 
         # Inverse of lengthscales should be drawn from truncated Gaussian | N(0, max_dist^2) |
@@ -59,26 +151,49 @@ class SpectralMixtureKernel(Kernel):
         # Draw means from Unif(0, 0.5 / minimum distance between two points)
         self.log_mixture_means.data.uniform_().mul_(0.5).div_(min_dist).log_()
         # Mixture weights should be roughly the stdv of the y values divided by the number of mixtures
-        self.log_mixture_weights.data.fill_(train_y.std() / self.n_mixtures).log_()
+        self.log_mixture_weights.data.fill_(train_y.std() / self.num_mixtures).log_()
 
     def forward(self, x1, x2):
-        batch_size, n, n_dims = x1.size()
+        batch_size, n, num_dims = x1.size()
         _, m, _ = x2.size()
-        if not n_dims == self.n_dims:
-            raise RuntimeError("The number of dimensions doesn't match what was supplied!")
 
-        mixture_weights = self.log_mixture_weights.view(self.n_mixtures, 1, 1, 1).exp()
-        mixture_means = self.log_mixture_means.view(self.n_mixtures, 1, 1, 1, self.n_dims).exp()
-        mixture_scales = self.log_mixture_scales.view(self.n_mixtures, 1, 1, 1, self.n_dims).exp()
-        distance = (x1.unsqueeze(-2) - x2.unsqueeze(-3)).abs()  # distance = x^(i) - z^(i)
+        if not num_dims == self.ard_num_dims:
+            raise RuntimeError(
+                "The SpectralMixtureKernel expected the input to have {} dimensionality "
+                "(based on the ard_num_dims argument). Got {}.".format(self.ard_num_dims, num_dims)
+            )
+        if not batch_size == self.batch_size:
+            raise RuntimeError(
+                "The SpectralMixtureKernel expected the input to have a batch_size of {} "
+                "(based on the batch_size argument). Got {}.".format(self.batch_size, batch_size)
+            )
+
+        # Expand x1 and x2 to account for the number of mixtures
+        # Should make x1/x2 (b x k x n x d) for k mixtures
+        x1_ = x1.unsqueeze(1)
+        x2_ = x2.unsqueeze(1)
+
+        # Compute distances - scaled by appropriate parameters
+        x1_exp = x1_ * self.mixture_scales
+        x2_exp = x2_ * self.mixture_scales
+        x1_cos = x1_ * self.mixture_means
+        x2_cos = x2_ * self.mixture_means
+
+        # Create grids
+        x1_exp_, x2_exp_ = self._create_input_grid(x1_exp, x2_exp)
+        x1_cos_, x2_cos_ = self._create_input_grid(x1_cos, x2_cos)
 
-        exp_term = (distance * mixture_scales).pow_(2).mul_(-2 * math.pi ** 2)
-        cos_term = (distance * mixture_means).mul_(2 * math.pi)
+        # Compute the exponential and cosine terms
+        exp_term = (x1_exp_ - x2_exp_).pow_(2).mul_(-2 * math.pi ** 2)
+        cos_term = (x1_cos_ - x2_cos_).mul_(2 * math.pi)
         res = exp_term.exp_() * cos_term.cos_()
 
-        # Product over dimensions
+        # Product omer dimensions
         res = res.prod(-1)
 
         # Sum over mixtures
-        res = (res * mixture_weights).sum(0)
+        mixture_weights = self.mixture_weights
+        while mixture_weights.dim() < res.dim():
+            mixture_weights.unsqueeze_(-1)
+        res = (res * mixture_weights).sum(1)
         return res
diff --git a/gpytorch/kernels/white_noise_kernel.py b/gpytorch/kernels/white_noise_kernel.py
index 27d98ea1..b89ac5b9 100644
--- a/gpytorch/kernels/white_noise_kernel.py
+++ b/gpytorch/kernels/white_noise_kernel.py
@@ -9,6 +9,30 @@ from gpytorch.lazy import DiagLazyTensor, ZeroLazyTensor
 
 
 class WhiteNoiseKernel(Kernel):
+    """
+    A "random" kernel that adds pre-specified white noise variance to training inputs.
+    This is most commonly used in conjunction with another kernel.
+
+    .. note::
+
+        The white noise is only applied to the portion of the kernel matrix that
+        represents the training data.
+
+    Args:
+        :attr:`variances` (Tensor `n x 1` or `b x n x 1`):
+            The random variances to be applied to training inputs.
+            `b` and `n` should correspond to the size of the training data.
+
+    Example:
+        >>> train_x = torch.randn(10, 5)
+        >>> wn_variances = torch.randn(10)
+        >>>
+        >>> covar_module = gpytorch.kernels.ScaleKernel(
+        >>>     gpytorch.kernels.WhiteNoiseKernel(wn_variances) + gpytorch.kernels.MaternKernel(nu=0.5)
+        >>> )
+        >>> covar = covar_module(train_x)  # Output: LazyVariable of size (10 x 10) (Matern kernel + random variances)
+    """
+
     def __init__(self, variances):
         super(WhiteNoiseKernel, self).__init__()
         self.register_buffer("variances", variances)


Test patch:
diff --git a/test/examples/test_kissgp_additive_regression.py b/test/examples/test_kissgp_additive_regression.py
index 6ae2b35e..864f623d 100644
--- a/test/examples/test_kissgp_additive_regression.py
+++ b/test/examples/test_kissgp_additive_regression.py
@@ -60,8 +60,8 @@ class TestKISSGPAdditiveRegression(unittest.TestCase):
             self.rng_state = torch.get_rng_state()
             torch.manual_seed(1)
             if torch.cuda.is_available():
-                torch.cuda.manual_seed_all(0)
-            random.seed(0)
+                torch.cuda.manual_seed_all(1)
+            random.seed(1)
 
     def tearDown(self):
         if hasattr(self, "rng_state"):
@@ -79,7 +79,7 @@ class TestKISSGPAdditiveRegression(unittest.TestCase):
 
             optimizer = optim.Adam(list(gp_model.parameters()) + list(likelihood.parameters()), lr=0.2)
             optimizer.n_iter = 0
-            for _ in range(20):
+            for _ in range(25):
                 optimizer.zero_grad()
                 output = gp_model(train_x)
                 loss = -mll(output, train_y)
diff --git a/test/examples/test_kronecker_multitask_gp_regression.py b/test/examples/test_kronecker_multitask_gp_regression.py
index 1f23567f..c16a6d87 100644
--- a/test/examples/test_kronecker_multitask_gp_regression.py
+++ b/test/examples/test_kronecker_multitask_gp_regression.py
@@ -33,7 +33,7 @@ class MultitaskGPModel(gpytorch.models.ExactGP):
         super(MultitaskGPModel, self).__init__(train_x, train_y, likelihood)
         self.mean_module = MultitaskMean(ConstantMean(), n_tasks=2)
         self.data_covar_module = RBFKernel()
-        self.covar_module = MultitaskKernel(self.data_covar_module, n_tasks=2, rank=1)
+        self.covar_module = MultitaskKernel(self.data_covar_module, n_tasks=2, rank=2)
 
     def forward(self, x):
         mean_x = self.mean_module(x)
diff --git a/test/examples/test_spectral_mixture_gp_regression.py b/test/examples/test_spectral_mixture_gp_regression.py
index a258d8d6..c2f521f8 100644
--- a/test/examples/test_spectral_mixture_gp_regression.py
+++ b/test/examples/test_spectral_mixture_gp_regression.py
@@ -29,11 +29,17 @@ test_y = torch.sin(test_x * (2 * pi))
 
 good_state_dict = OrderedDict(
     [
-        ("likelihood.log_noise", torch.tensor([-13.4054])),
+        ("likelihood.log_noise", torch.tensor([-5])),
         ("mean_module.constant", torch.tensor([[0.4615]])),
-        ("covar_module.log_mixture_weights", torch.tensor([-0.7277, -15.1212, -0.5511, -6.3787])),
-        ("covar_module.log_mixture_means", torch.tensor([[-0.1201], [0.6013], [-3.7319], [0.2380]])),
-        ("covar_module.log_mixture_scales", torch.tensor([[-1.9713], [2.6217], [-3.9268], [-4.7071]])),
+        ("covar_module.log_mixture_weights", torch.tensor([-0.7277, -15.1212, -0.5511, -6.3787]).unsqueeze(0)),
+        (
+            "covar_module.log_mixture_means",
+            torch.tensor([[-0.1201], [0.6013], [-3.7319], [0.2380]]).unsqueeze(0).unsqueeze(-2),
+        ),
+        (
+            "covar_module.log_mixture_scales",
+            torch.tensor([[-1.9713], [2.6217], [-3.9268], [-4.7071]]).unsqueeze(0).unsqueeze(-2),
+        ),
     ]
 )
 
@@ -42,7 +48,7 @@ class SpectralMixtureGPModel(gpytorch.models.ExactGP):
     def __init__(self, train_x, train_y, likelihood):
         super(SpectralMixtureGPModel, self).__init__(train_x, train_y, likelihood)
         self.mean_module = ConstantMean(prior=SmoothedBoxPrior(-1, 1))
-        self.covar_module = SpectralMixtureKernel(n_mixtures=4)
+        self.covar_module = SpectralMixtureKernel(num_mixtures=4, ard_num_dims=1)
         self.covar_module.initialize_from_data(train_x, train_y)
 
     def forward(self, x):
@@ -54,11 +60,12 @@ class SpectralMixtureGPModel(gpytorch.models.ExactGP):
 class TestSpectralMixtureGPRegression(unittest.TestCase):
     def setUp(self):
         if os.getenv("UNLOCK_SEED") is None or os.getenv("UNLOCK_SEED").lower() == "false":
+            seed = 4
             self.rng_state = torch.get_rng_state()
-            torch.manual_seed(4)
+            torch.manual_seed(seed)
             if torch.cuda.is_available():
-                torch.cuda.manual_seed_all(4)
-            random.seed(4)
+                torch.cuda.manual_seed_all(seed)
+            random.seed(seed)
 
     def tearDown(self):
         if hasattr(self, "rng_state"):
@@ -97,6 +104,7 @@ class TestSpectralMixtureGPRegression(unittest.TestCase):
             gp_model.load_state_dict(good_state_dict, strict=False)
 
             # Test the model
+        with torch.no_grad(), gpytorch.settings.max_cg_iterations(100):
             gp_model.eval()
             likelihood.eval()
             test_preds = likelihood(gp_model(test_x)).mean()
diff --git a/test/kernels/test_cosine_kernel.py b/test/kernels/test_cosine_kernel.py
new file mode 100644
index 00000000..e76f6e64
--- /dev/null
+++ b/test/kernels/test_cosine_kernel.py
@@ -0,0 +1,31 @@
+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
+from torch.autograd import Variable
+from gpytorch.kernels import CosineKernel
+
+
+class TestCosineKernel(unittest.TestCase):
+    def test_computes_periodic_function(self):
+        a = torch.Tensor([4, 2, 8]).view(3, 1)
+        b = torch.Tensor([0, 2]).view(2, 1)
+        period = 1
+        kernel = CosineKernel().initialize(log_period_length=math.log(period))
+        kernel.eval()
+
+        actual = torch.zeros(3, 2)
+        for i in range(3):
+            for j in range(2):
+                actual[i, j] = torch.cos(math.pi * (a[i] - b[j]) / period)
+
+        res = kernel(Variable(a), Variable(b)).evaluate()
+        self.assertLess(torch.norm(res - actual), 1e-5)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/kernels/test_matern_kernel.py b/test/kernels/test_matern_kernel.py
new file mode 100644
index 00000000..2753bc22
--- /dev/null
+++ b/test/kernels/test_matern_kernel.py
@@ -0,0 +1,53 @@
+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
+from gpytorch.kernels import MaternKernel
+
+
+class TestMaternKernel(unittest.TestCase):
+    def test_forward_nu_1_over_2(self):
+        a = torch.Tensor([4, 2, 8]).view(3, 1)
+        b = torch.Tensor([0, 2]).view(2, 1)
+        lengthscale = 2
+
+        kernel = MaternKernel(nu=0.5).initialize(log_lengthscale=math.log(lengthscale))
+        kernel.eval()
+
+        actual = torch.Tensor([[4, 2], [2, 0], [8, 6]]).div_(-lengthscale).exp()
+        res = kernel(a, b).evaluate()
+        self.assertLess(torch.norm(res - actual), 1e-3)
+
+    def test_forward_nu_3_over_2(self):
+        a = torch.Tensor([4, 2, 8]).view(3, 1)
+        b = torch.Tensor([0, 2]).view(2, 1)
+        lengthscale = 2
+
+        kernel = MaternKernel(nu=1.5).initialize(log_lengthscale=math.log(lengthscale))
+        kernel.eval()
+
+        dist = torch.Tensor([[4, 2], [2, 0], [8, 6]]).mul_(math.sqrt(3) / lengthscale)
+        actual = (dist + 1).mul(torch.exp(-dist))
+        res = kernel(a, b).evaluate()
+        self.assertLess(torch.norm(res - actual), 1e-3)
+
+    def test_forward_nu_5_over_2(self):
+        a = torch.Tensor([4, 2, 8]).view(3, 1)
+        b = torch.Tensor([0, 2]).view(2, 1)
+        lengthscale = 2
+
+        kernel = MaternKernel(nu=2.5).initialize(log_lengthscale=math.log(lengthscale))
+        kernel.eval()
+
+        dist = torch.Tensor([[4, 2], [2, 0], [8, 6]]).mul_(math.sqrt(5) / lengthscale)
+        actual = (dist ** 2 / 3 + dist + 1).mul(torch.exp(-dist))
+        res = kernel(a, b).evaluate()
+        self.assertLess(torch.norm(res - actual), 1e-3)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/kernels/test_periodic_kernel.py b/test/kernels/test_periodic_kernel.py
index f3719081..88b2d7d8 100644
--- a/test/kernels/test_periodic_kernel.py
+++ b/test/kernels/test_periodic_kernel.py
@@ -14,17 +14,17 @@ class TestPeriodicKernel(unittest.TestCase):
         a = torch.Tensor([4, 2, 8]).view(3, 1)
         b = torch.Tensor([0, 2]).view(2, 1)
         lengthscale = 2
-        period = 1
+        period = 3
         kernel = PeriodicKernel().initialize(log_lengthscale=math.log(lengthscale), log_period_length=math.log(period))
         kernel.eval()
 
         actual = torch.zeros(3, 2)
         for i in range(3):
             for j in range(2):
-                val = 2 * torch.pow(torch.sin(math.pi * (a[i] - b[j])), 2) / lengthscale
-                actual[i, j] = torch.exp(-val)[0]
+                val = 2 * torch.pow(torch.sin(math.pi * (a[i] - b[j]) / 3), 2) / lengthscale
+                actual[i, j] = torch.exp(-val).item()
 
-        res = kernel(a, b).evaluate().data
+        res = kernel(a, b).evaluate()
         self.assertLess(torch.norm(res - actual), 1e-5)
 
 
diff --git a/test/kernels/test_spectral_mixture_kernel.py b/test/kernels/test_spectral_mixture_kernel.py
new file mode 100644
index 00000000..dc1f68e0
--- /dev/null
+++ b/test/kernels/test_spectral_mixture_kernel.py
@@ -0,0 +1,42 @@
+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
+from torch.autograd import Variable
+from gpytorch.kernels import SpectralMixtureKernel
+
+
+class TestSpectralMixtureKernel(unittest.TestCase):
+    def test_computes_periodic_function(self):
+        a = torch.Tensor([4, 2, 8]).view(3, 1)
+        b = torch.Tensor([0, 2]).view(2, 1)
+        means = [1, 2]
+        scales = [0.5, 0.25]
+        weights = [4, 2]
+        kernel = SpectralMixtureKernel(num_mixtures=2)
+        kernel.initialize(
+            log_mixture_weights=torch.Tensor([[4, 2]]).log(),
+            log_mixture_means=torch.Tensor([[[[1]], [[2]]]]).log(),
+            log_mixture_scales=torch.Tensor([[[[0.5]], [[0.25]]]]).log(),
+        )
+        kernel.eval()
+
+        actual = torch.zeros(3, 2)
+        for i in range(3):
+            for j in range(2):
+                for k in range(2):
+                    new_term = torch.cos(2 * math.pi * (a[i] - b[j]) * means[k])
+                    new_term *= torch.exp(-2 * (math.pi * (a[i] - b[j])) ** 2 * scales[k] ** 2)
+                    new_term *= weights[k]
+                    actual[i, j] += new_term.item()
+
+        res = kernel(Variable(a), Variable(b)).evaluate()
+        self.assertLess(torch.norm(res - actual), 1e-5)
+
+
+if __name__ == "__main__":
+    unittest.main()