Repository: cornellius-gp/gpytorch
Pull request #450 (2018-12-31 00:32:54)

Support for multi-device GPs
This pull request consists of a new `CatLazyTensor` and a new `DistributedKernel`. Together they allow for GP training and inference on multiple devices. 

`CatLazyTensor` allows the user to concatenate `LazyTensors` with each possibly on a different device. It then distributes `matmul` calls across each device as necessary. `CatLazyTensor` also can be sliced along any dimension including batch dimensions. This means one can concatenate multiple types of `LazyTensors` together and obtain a heterogeneous `LazyTensor` containing slices from the underlying `LazyTensors`.

`DistributedKernel` handles creating distributed `CatLazyTensor`s in its `forward` method. 

Included some notebooks of some tests used during development and a Distributed GP Regression notebook example. 

TODOs: The test notebooks should be removed and turned into unittests. Also some more tests should be written.

jacobrgardner: NB: the for loops in `CatLazyTensor._matmul` and elsewhere are correct: PyTorch does these calls asynchronously (see https://discuss.pytorch.org/t/matmul-on-multiple-gpus/33122/3).
KeAWang: We might also want to do a `gpytorch.cat` function much like a `torch.cat` function
jacobrgardner: @KeAWang @gpleiss I rewrote `MultiDeviceKernel` to inherit from `torch.nn.parallel.DataParallel`. We now call the kernels in parallel, and rely on the upstream code for splitting `x1` and distributing the kernel module on multiple devices.

MultiDeviceKernel can't incorporate the CPU at all now as a result (as `DataParallel` only supports multiple GPUs). However, this makes the multi GPU setting more efficient, and I can also confirm that this all does indeed work with multiple GPUs!

I've also updated the example notebook.

Patch:
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 41c57de9..fb4e62f4 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -41,6 +41,10 @@ class _SubDistribution(object):
     pass
 
 
+class _Kernel(object):
+    pass
+
+
 # Putting all of our dirty hacks together
 class Mock(MagicMock):
     __metaclass__ = type
@@ -53,6 +57,8 @@ class Mock(MagicMock):
             return _Distribution
         elif "Normal" in name or "Gamma" in name or "Wishart" in name:
             return _SubDistribution
+        elif "Kernel" in name or "Parallel" in name:
+            return _Kernel
         else:
             res = MagicMock()
             res.Module = ModuleMock
@@ -68,6 +74,7 @@ MOCK_MODULES = [
     "torch.autograd",
     "torch.nn",
     "torch.nn.functional",
+    "torch.nn.parallel",
     "torch.optim",
     "torch.utils",
     "torch.utils.data",
diff --git a/docs/source/kernels.rst b/docs/source/kernels.rst
index 0d4375a7..5ac842d3 100644
--- a/docs/source/kernels.rst
+++ b/docs/source/kernels.rst
@@ -76,6 +76,13 @@ Composition/Decoration Kernels
 :hidden:`AdditiveStructureKernel`
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+:hidden:`MultiDeviceKernel`
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. autoclass:: MultiDeviceKernel
+   :members:
+
+
 .. autoclass:: AdditiveStructureKernel
    :members:
 
diff --git a/docs/source/lazy.rst b/docs/source/lazy.rst
index 733d2afa..e1167850 100644
--- a/docs/source/lazy.rst
+++ b/docs/source/lazy.rst
@@ -34,6 +34,13 @@ Structured LazyTensors
 .. autoclass:: BlockDiagLazyTensor
    :members:
 
+:hidden:`CatLazyTensor`
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. autoclass:: CatLazyTensor
+   :members:
+
+
 :hidden:`CholLazyTensor`
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/examples/01_Simple_GP_Regression/Simple_GP_Regression.ipynb b/examples/01_Simple_GP_Regression/Simple_GP_Regression.ipynb
index da7b16c1..aab5fcef 100644
--- a/examples/01_Simple_GP_Regression/Simple_GP_Regression.ipynb
+++ b/examples/01_Simple_GP_Regression/Simple_GP_Regression.ipynb
@@ -358,7 +358,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.0"
+   "version": "3.6.7"
   }
  },
  "nbformat": 4,
diff --git a/examples/01_Simple_GP_Regression/Simple_MultiGPU_GP_Regression.ipynb b/examples/01_Simple_GP_Regression/Simple_MultiGPU_GP_Regression.ipynb
new file mode 100644
index 00000000..5d67b429
--- /dev/null
+++ b/examples/01_Simple_GP_Regression/Simple_MultiGPU_GP_Regression.ipynb
@@ -0,0 +1,260 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Exact GP Regression with Multiple GPUs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "import math\n",
+    "import torch\n",
+    "import gpytorch\n",
+    "from matplotlib import pyplot as plt\n",
+    "\n",
+    "%matplotlib inline\n",
+    "%load_ext autoreload\n",
+    "%autoreload 2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Training data is 11 points in [0,1] inclusive regularly spaced\n",
+    "train_x = torch.linspace(0, 1, 100)\n",
+    "# True function is sin(2*pi*x) with Gaussian noise\n",
+    "train_y = torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * 0.2\n",
+    "\n",
+    "train_x = train_x.cuda()\n",
+    "train_y = train_y.cuda()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# We will use the simplest form of GP model, exact inference\n",
+    "class ExactGPModel(gpytorch.models.ExactGP):\n",
+    "    def __init__(self, train_x, train_y, likelihood):\n",
+    "        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)\n",
+    "        self.mean_module = gpytorch.means.ConstantMean()\n",
+    "        base_covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\n",
+    "        \n",
+    "        devices = [0, 1]\n",
+    "        \n",
+    "        self.covar_module = gpytorch.kernels.MultiDeviceKernel(base_covar_module, device_ids=devices, output_device=torch.device('cuda', 0))\n",
+    "    \n",
+    "    def forward(self, x):\n",
+    "        mean_x = self.mean_module(x)\n",
+    "        covar_x = self.covar_module(x)\n",
+    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
+    "\n",
+    "likelihood = gpytorch.likelihoods.GaussianLikelihood().to(\"cuda:0\")\n",
+    "model = ExactGPModel(train_x, train_y, likelihood).to(\"cuda:0\")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/jrg365/gpytorch/gpytorch/utils/pivoted_cholesky.py:103: UserWarning: torch.potrs is deprecated in favour of torch.cholesky_solve and will be removed in the next release. Please use torch.cholesky instead and note that the :attr:`upper` argument in torch.cholesky_solve defaults to ``False``.\n",
+      "  R = torch.potrs(low_rank_mat, torch.cholesky(shifted_mat, upper=True))\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Iter 1/50 - Loss: 0.937   log_lengthscale: -0.367   log_noise: -0.367\n",
+      "Iter 2/50 - Loss: 0.906   log_lengthscale: -0.439   log_noise: -0.439\n",
+      "Iter 3/50 - Loss: 0.873   log_lengthscale: -0.513   log_noise: -0.514\n",
+      "Iter 4/50 - Loss: 0.836   log_lengthscale: -0.589   log_noise: -0.590\n",
+      "Iter 5/50 - Loss: 0.794   log_lengthscale: -0.666   log_noise: -0.668\n",
+      "Iter 6/50 - Loss: 0.748   log_lengthscale: -0.744   log_noise: -0.747\n",
+      "Iter 7/50 - Loss: 0.699   log_lengthscale: -0.823   log_noise: -0.828\n",
+      "Iter 8/50 - Loss: 0.649   log_lengthscale: -0.905   log_noise: -0.911\n",
+      "Iter 9/50 - Loss: 0.601   log_lengthscale: -0.989   log_noise: -0.996\n",
+      "Iter 10/50 - Loss: 0.557   log_lengthscale: -1.072   log_noise: -1.082\n",
+      "Iter 11/50 - Loss: 0.518   log_lengthscale: -1.153   log_noise: -1.169\n",
+      "Iter 12/50 - Loss: 0.483   log_lengthscale: -1.230   log_noise: -1.258\n",
+      "Iter 13/50 - Loss: 0.447   log_lengthscale: -1.300   log_noise: -1.349\n",
+      "Iter 14/50 - Loss: 0.413   log_lengthscale: -1.360   log_noise: -1.440\n",
+      "Iter 15/50 - Loss: 0.387   log_lengthscale: -1.412   log_noise: -1.532\n",
+      "Iter 16/50 - Loss: 0.350   log_lengthscale: -1.449   log_noise: -1.625\n",
+      "Iter 17/50 - Loss: 0.321   log_lengthscale: -1.477   log_noise: -1.719\n",
+      "Iter 18/50 - Loss: 0.297   log_lengthscale: -1.499   log_noise: -1.813\n",
+      "Iter 19/50 - Loss: 0.257   log_lengthscale: -1.510   log_noise: -1.908\n",
+      "Iter 20/50 - Loss: 0.243   log_lengthscale: -1.513   log_noise: -2.002\n",
+      "Iter 21/50 - Loss: 0.208   log_lengthscale: -1.506   log_noise: -2.096\n",
+      "Iter 22/50 - Loss: 0.178   log_lengthscale: -1.491   log_noise: -2.190\n",
+      "Iter 23/50 - Loss: 0.154   log_lengthscale: -1.468   log_noise: -2.283\n",
+      "Iter 24/50 - Loss: 0.137   log_lengthscale: -1.444   log_noise: -2.375\n",
+      "Iter 25/50 - Loss: 0.113   log_lengthscale: -1.413   log_noise: -2.465\n",
+      "Iter 26/50 - Loss: 0.088   log_lengthscale: -1.376   log_noise: -2.554\n",
+      "Iter 27/50 - Loss: 0.074   log_lengthscale: -1.338   log_noise: -2.640\n",
+      "Iter 28/50 - Loss: 0.068   log_lengthscale: -1.302   log_noise: -2.723\n",
+      "Iter 29/50 - Loss: 0.054   log_lengthscale: -1.263   log_noise: -2.803\n",
+      "Iter 30/50 - Loss: 0.047   log_lengthscale: -1.225   log_noise: -2.879\n",
+      "Iter 31/50 - Loss: 0.043   log_lengthscale: -1.190   log_noise: -2.950\n",
+      "Iter 32/50 - Loss: 0.045   log_lengthscale: -1.162   log_noise: -3.016\n",
+      "Iter 33/50 - Loss: 0.048   log_lengthscale: -1.139   log_noise: -3.077\n",
+      "Iter 34/50 - Loss: 0.052   log_lengthscale: -1.124   log_noise: -3.130\n",
+      "Iter 35/50 - Loss: 0.057   log_lengthscale: -1.117   log_noise: -3.177\n",
+      "Iter 36/50 - Loss: 0.060   log_lengthscale: -1.114   log_noise: -3.216\n",
+      "Iter 37/50 - Loss: 0.062   log_lengthscale: -1.121   log_noise: -3.247\n",
+      "Iter 38/50 - Loss: 0.065   log_lengthscale: -1.135   log_noise: -3.270\n",
+      "Iter 39/50 - Loss: 0.068   log_lengthscale: -1.155   log_noise: -3.285\n",
+      "Iter 40/50 - Loss: 0.065   log_lengthscale: -1.182   log_noise: -3.293\n",
+      "Iter 41/50 - Loss: 0.066   log_lengthscale: -1.210   log_noise: -3.294\n",
+      "Iter 42/50 - Loss: 0.063   log_lengthscale: -1.237   log_noise: -3.289\n",
+      "Iter 43/50 - Loss: 0.061   log_lengthscale: -1.264   log_noise: -3.277\n",
+      "Iter 44/50 - Loss: 0.058   log_lengthscale: -1.290   log_noise: -3.261\n",
+      "Iter 45/50 - Loss: 0.057   log_lengthscale: -1.309   log_noise: -3.240\n",
+      "Iter 46/50 - Loss: 0.054   log_lengthscale: -1.326   log_noise: -3.216\n",
+      "Iter 47/50 - Loss: 0.057   log_lengthscale: -1.340   log_noise: -3.190\n",
+      "Iter 48/50 - Loss: 0.051   log_lengthscale: -1.347   log_noise: -3.162\n",
+      "Iter 49/50 - Loss: 0.052   log_lengthscale: -1.351   log_noise: -3.133\n",
+      "Iter 50/50 - Loss: 0.047   log_lengthscale: -1.349   log_noise: -3.104\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Find optimal model hyperparameters\n",
+    "model.train()\n",
+    "likelihood.train()\n",
+    "\n",
+    "# Use the adam optimizer\n",
+    "optimizer = torch.optim.Adam([\n",
+    "    {'params': model.parameters()},  # Includes GaussianLikelihood parameters\n",
+    "], lr=0.1)\n",
+    "\n",
+    "# \"Loss\" for GPs - the marginal log likelihood\n",
+    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
+    "\n",
+    "with gpytorch.settings.max_preconditioner_size(5):\n",
+    "    training_iter = 50\n",
+    "    for i in range(training_iter):\n",
+    "        # Zero gradients from previous iteration\n",
+    "        optimizer.zero_grad()\n",
+    "        # Output from model\n",
+    "        output = model(train_x)\n",
+    "        # Calc loss and backprop gradients\n",
+    "        loss = -mll(output, train_y)\n",
+    "        loss.backward()\n",
+    "        print('Iter %d/%d - Loss: %.3f   log_lengthscale: %.3f   log_noise: %.3f' % (\n",
+    "            i + 1, training_iter, loss.item(),\n",
+    "            model.covar_module.module.base_kernel.log_lengthscale.item(),\n",
+    "            model.likelihood.log_noise.item()\n",
+    "        ))\n",
+    "        optimizer.step()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/jrg365/gpytorch/gpytorch/utils/pivoted_cholesky.py:103: UserWarning: torch.potrs is deprecated in favour of torch.cholesky_solve and will be removed in the next release. Please use torch.cholesky instead and note that the :attr:`upper` argument in torch.cholesky_solve defaults to ``False``.\n",
+      "  R = torch.potrs(low_rank_mat, torch.cholesky(shifted_mat, upper=True))\n"
+     ]
+    }
+   ],
+   "source": [
+    "# Get into evaluation (predictive posterior) mode\n",
+    "model.eval()\n",
+    "likelihood.eval()\n",
+    "\n",
+    "# Test points are regularly spaced along [0,1]\n",
+    "# Make predictions by feeding model through likelihood\n",
+    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
+    "    test_x = torch.linspace(0, 1, 51)\n",
+    "    observed_pred = model(test_x)\n",
+    "    observed_pred = likelihood(observed_pred)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAADGCAYAAAAwqi48AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJztnXl4U1X6+D8nadp0pQstLZS1LLK0ZRcEQUBgREBBURYdBR0FBsWf4qgjI7jNjBsjCjqDjKAOgoKgfgEXUJS1skjZpUApUOi+b2na5vz+SBPbJmlTGkoK5/M8fZrcnHvue29y33vOe95FSClRKBSKqmiutgAKhcL9UIpBoVDYoBSDQqGwQSkGhUJhg1IMCoXCBqUYFAqFDQ1WDEIIvRBirxDikBDimBDiRVcIplAorh6ioX4MQggB+EopC4UQOmAnMFdKGecKARUKRePj0dAOpFmzFFa+1VX+Ka8phaIJ4xIbgxBCK4SIB9KBLVLKX1zRr0KhuDo0eMQAIKWsAHoKIQKBDUKIHlLKo1XbCCEeAR4B8PX17XPDDTe44tAKhaIeHDhwIFNKGVpXuwbbGGw6FGIBUCSlfNNRm759+8r9+/e79LgKhaJuhBAHpJR962rnilWJ0MqRAkIIb+BW4LeG9qtQKK4erphKRAAfCSG0mBXN51LKjS7oV6FQXCVcsSpxGOjlAlkUCoWb4BLjo+LapaysjOTkZAwGw9UWRVEP9Ho9kZGR6HS6y9pfKQZFrSQnJ+Pv70+7du0w+7Ip3B0pJVlZWSQnJ9O+ffvL6kPFSihqxWAwEBISopRCE0IIQUhISINGeUoxKOpEKYWmR0O/M6UYFG5PcnIyd9xxB506dSIqKoq5c+diNBoBWLlyJXPmzLnKEtri5+dnd7tWq6Vnz550796d2NhYFi1ahMlkqrWvpKQkPv300yshpkOUYlC4nJSUFIYOHUpqamqD+5JSMnHiRO68805OnTpFQkIChYWFPP/88y6Q1D7l5eVXrG9vb2/i4+M5duwYW7ZsYfPmzbz4Yu0ByVdDMSClbPS/Pn36SEXT4Pjx4/XeZ9asWVKj0chZs2Y1+Phbt26VN998c7VteXl5Mjg4WBYVFckVK1bI8ePHy9GjR8vOnTvLhQsXSimlLCwslGPGjJExMTGye/fucs2aNVJKKffv3y+HDBkie/fuLUeNGiUvXbokpZRy6NCh8rnnnpNDhgyRCxculG3btpUVFRVSSimLiopkZGSkNBqN8vTp03L06NGyd+/ecvDgwfLEiRNSSikTExPlgAEDZN++feX8+fOlr6+v3fOpuf3MmTMyODhYmkwmefbsWTl48GDZq1cv2atXL7lr1y4ppZQ33nijDAgIkLGxsXLRokUO29XE3ncH7JdO3KNKMShqpT6KQa/XS8yRtdX+9Hr9ZR9/8eLF8oknnrDZ3rNnT3no0CG5YsUKGR4eLjMzM2VxcbHs3r273Ldvn1y3bp18+OGHre1zc3Ol0WiUAwcOlOnp6VJKKdesWSOnT58upTQrhqqKbPz48fLHH3+0tnvooYeklFIOHz5cJiQkSCmljIuLk8OGDZNSSjlu3Dj50UcfSSmlXLJkidOKQUopAwMDZWpqqiwqKpIlJSVSSikTEhKk5T7Ztm2bvP32263tHbWrSUMUg5pKKFxGYmIiU6dOxcfHBwAfHx+mTZvG2bNnL7tPKaVdQ1rV7SNHjiQkJARvb28mTpzIzp07iY6OZuvWrTzzzDPs2LGDZs2acfLkSY4ePcrIkSPp2bMnr7zyCsnJydY+77333mqvP/vsMwDWrFnDvffeS2FhIbt372bSpEn07NmTRx99lJSUFAB27drFlClTALj//vvrfY5g9hn505/+RHR0NJMmTeL48eN22zvbriEoPwaFy4iIiCAgIACDwYBer8dgMBAQEEB4ePhl99m9e3e++OKLatvy8/O5cOECUVFRHDhwwEZxCCHo3LkzBw4cYPPmzTz33HOMGjWKCRMm0L17d/bs2WP3WL6+vtbX48eP57nnniM7O5sDBw4wfPhwioqKCAwMJD4+3u7+l7MSkJiYiFarJSwsjBdffJEWLVpw6NAhTCYTer3e7j7/+te/nGrXENSIQeFS0tLSmDlzJnFxccycObPBBsgRI0ZQXFzMxx9/DEBFRQVPPfUUDz74oHVksmXLFrKzsykpKeHLL79k0KBBXLp0CR8fH+677z7mzZvHr7/+SpcuXcjIyLAqhrKyMo4dO2b3uH5+fvTv35+5c+cyduxYtFotAQEBtG/fnrVr1wLmJ/2hQ4cAGDRoEGvWrAFg1apVTp1bRkYGM2fOZM6cOQghyMvLIyIiAo1GwyeffEJFRQUA/v7+FBQUWPdz1M6lODPfcPWfsjE0HS7H+Ohqzp8/L8eOHSs7duwoO3ToIOfMmSMNBoOUUsoVK1bISZMmyTFjxlQzPn777bcyOjpaxsbGyr59+8p9+/ZJKaU8ePCgvPnmm2VMTIzs1q2bXLZsmZTSbGOwtLGwdu1aCciffvrJui0xMVGOHj1axsTEyK5du8oXX3zRut1ifPzHP/7h0Mag0WhkbGys7Natm4yJiZFvvPGG1ciZkJAgo6Oj5Y033iifffZZax9Go1EOHz5cxsTEyEWLFjlsV5OG2Bhcno/BGVQ+hqbDiRMn6Nq169UWQ3EZ2PvuGi0fg0KhuPZQikGhUNigFINCobBBKQaFQmGDUgwKhcIGVySDbS2E2CaEOFFZom6uKwRTKBRXD1eMGMqBp6SUXYEBwJ+FEN1c0K9CAZg9Cqu6GZeXlxMaGsrYsWOvolTXNg1WDFLKFCnlr5WvC4ATQKuG9qtQWPD19eXo0aOUlJQAZk/HVq3UT+xK4lIbgxCiHeaM0apEncKl3HbbbWzatAmA1atXWwOWAIqKipgxYwb9+vWjV69efPXVV4A5j8HNN99M79696d27N7t37wbgp59+4pZbbuHuu+/mhhtuYNq0aVwNRz93xmVBVEIIP+AL4AkpZb6dz60l6tq0aeOqwyoakSeeAAfxQ5dNz57w9tt1t5s8eTIvvfQSY8eO5fDhw8yYMYMdO3YA8OqrrzJ8+HA+/PBDcnNz6d+/P7feeithYWFs2bIFvV7PqVOnmDJlChaP24MHD3Ls2DFatmzJoEGD2LVrF4MHD3btyTVhXKIYhBA6zEphlZRyvb02UsplwDIwu0S74riK64eYmBiSkpJYvXo1Y8aMqfbZ999/z9dff82bb5qrIhoMBs6fP0/Lli2ZM2cO8fHxaLVaEhISrPv079+fyMhIAHr27ElSUpJSDFVosGIQ5ljT/wInpJSLGi6Swl1x5sl+JRk/fjzz5s3jp59+Iisry7pdSskXX3xBly5dqrVfuHChw/BkLy8v62utVntF07k1RVxhYxgE3A8MF0LEV/6NqWsnhaK+zJgxgxdeeIHo6Ohq20ePHs27775rtRMcPHgQaKTw5GsUV6xK7JRSCilljJSyZ+XfZlcIp1BUJTIykrlzbd1k/va3v1FWVkZMTAw9evTgb3/7GwCzZ8/mo48+YsCAASQkJFRLxKKoHRV2ragVFXbddFFh1wqFwqUoxaBQKGxQikGhUNigFINCobBBKQaFQmGDUgwKhcIGpRgUTYLU1FQmT55MVFQU3bp1Y8yYMdVcnJ1lx44ddO/enZ49e3Lx4kXuvvtuu+1uueUWrucldVWJSlEv/rWl/jdjbfy/kZ3rbCOlZMKECTzwwAPWoi7x8fGkpaXRuXPd+1dl1apVzJs3j+nTpwOwbt26+gt9HaBGDC5GSkmxsZyMglKSMotIzCgkKbOI81nFXMguJrvISHmF6WqL2aTYtm0bOp2OmTNnWrf17NmTwYMH8/TTT9OjRw+io6OttSYdhVUvX76czz//nJdeeolp06aRlJREjx49ACgpKWHy5MnExMRw7733WnM/gDlIa+DAgfTu3ZtJkyZRWFgIQLt27ViwYAG9e/cmOjqa3377DYDCwkKmT59OdHQ0MTEx1hJ7jvpxR9SIoQFUmCTpBQZS8wyk5Zv/5xvKqTDV7k0qBPh6etDMR0eYvxeRQT5EBnmj12kbSfKmxdGjR+nTp4/N9vXr1xMfH8+hQ4fIzMykX79+DBkyBLAfVv3www+zc+dOxo4dy913301SUpK1r/fffx8fHx8OHz7M4cOH6d27NwCZmZm88sorbN26FV9fX1577TUWLVrECy+8AEDz5s359ddfee+993jzzTdZvnw5L7/8Ms2aNePIkSMA5OTk1NmPu6EUQz0pLa/gbGYRp9MLOZdVjLH896d/flY6H//9Sf74/L8ICA512IeUUFhazqVLl3iusn2zkFBC/LyIau7LDREBBPt6NsbpNGl27tzJlClT0Gq1tGjRgqFDh7Jv3z4CAgLqHVa9fft2Hn/8ccAc4h0TEwNAXFwcx48fZ9CgQQAYjUYGDhxo3W/ixIkA9OnTh/XrzRkHtm7dap3yAAQFBbFx48Za+3E3lGJwAikl57KKOZScy7msYocjgu9XvcfZo/v5/n9LufvxhXX2W7N9ZkEpmQWl/HI2mxYBerqE+9MtIgBvz+t7JNG9e3e7toDa4nwuJ6zaXrVqKSUjR45k9erVtR6n6jGklDZ91dWPu6FsDLVgKKvgwLlsVu5OYsPBiyRmFNlVCn8ZG8OTo7qwe+NqpJTs3riaJ0d14S9jY+z2W1v7/Kx0ljx1H6eSLrA9IYMPd51l28l08krKrvTpui3Dhw+ntLSUDz74wLpt3759BAUF8dlnn1FRUUFGRgbbt2+nf//+l3WMIUOGWKtUHz16lMOHDwMwYMAAdu3axenTpwEoLi6uczVk1KhRLFmyxPo+Jyfnsvq5mijFYIcSYwW7Tmfy351n2Z6QSW5x7Tfl/I+20nvYWHRe5kQgHp5e+DULZu47nzvVXuelJ3rQSFp1uIGNH75lHUUAGMtNxJ/PZeWuJL45kkJ2kdGFZ9o0EEKwYcMGtmzZQlRUFN27d2fhwoVMnTqVmJgYYmNjGT58OK+//jrh4eGXdYxZs2ZRWFhITEwMr7/+ulXBhIaGsnLlSqZMmUJMTAwDBgywGhkdMX/+fHJycujRowexsbFs27btsvq5mqiw6yqYRwg5/ByfwH9feqJOW0FV1i1ewJ7Nn6HVeVJuLAXgprFTHE4pqravKDMiwWx8qIGHpxevbzxsfa8RgujIAAZ0CMHH88rPBFXYddNFhV03kAqT5NfzOXy46yx7z2az6aMl1Z7azlCQmwVCWJUCYDOlsEwT8rMzKMjN4qaxU5Amk3mubEcp9B4+jvkf/1Btm0lKDl3IY8WuJPYlZaulT8UV4bo3PiZlFvFzQgbZRUb+MjbG5sbevXG1zVPbHtMXLCE/K52vl73Gkd1bKSs14OHphd7bl2l/fYslT91HcHgrq8KZvsA8Bx05dVa1fQCERgNSovfxczhiMZab2Hkqk+OX8hnRNYzIIB8XXRGF4joZMaSkpDB06FBSU1Ot2/INZXwVf5ENBy9a5+1z316Db7NgPDzNlmadl97uU9sRASFh6H38KDeW4uHpRbmxlMK8bP7z7AwSj+xj/5YvbYyNVfcRwvx1xAweRd9b7yR++zfkZ2fUeszsIiPrDiTzw4k0SstVTkOFa3CJYhBCfCiESBdCHHVFf67m5ZdfZufOnbz00ktIaZ42fLLnHIkZRdXa7dn8OUV52dVubL2PH0hpnQLUhb0pRV3TBMu04sn3NjBo3FRMFRXovPSUFOQ5NZ2REg4n5/HJnnOczSyqs319UcVYmh4N/c5cYnwUQgwBCoGPpZQ96mrfWMZHb29vDAaDzfaaU4OaUwgLQqOhx8AR+Ac1Z8+mNQy8fbJdY2J+VjofvvgKg8cvorQ4hPRkI8fiTpKTZkJKLebynhWV/zOAM3j7ZjHt2dlExfjg5f37d+BIFmemM2D2qoxtHcjNHZvjoW243j979iz+/v6EhITYXedXuB9SSrKysigoKKB9+/bVPnPW+OiyVYnK8nQb3UkxpKSkMG/ePL788kuKi4uty4LjH3mm2ty9pm3A0u7wju8oL7NdHvTw9OKZD45zLM6PpOPenNhroLTk9/40GomHZzZGwxmgFLMpxwPQ4eHZhnJjcLW2kZ0MRMWUEBVTQmjkOb77+B82stSUuS6a+3txW49wmvt51d24FsrKykhOTrarYBXui16vJzIyEp1OV227s4qh0YyPV6NEXUREBJ7evpSUlFSbGtS8wezZBvQ+fsxd/Bn/fu4hSkuKzJ/pYgmNnIeUd/LqA36VeycBezGX69wLnEJo8+nSZzABwaGkXUikMCcLv6AQzhzeS7nRBPgA7YAoTKb+nE8YxsUzN7FtbTA6rwiahTxLWakJre5nyo3FtRohHZFZUMrqX84ztEsoMZGBl30NdTqdzVNHce3TaIqhsUvUSSk5cC6HX08mcdPYKQwYcy9xmz9zaCewzPOrtjPbHAqAiSAeo7xsEClnoXXnfILD32bC7O4c3PaO00/36iOT4+i8EokeBOMfGYbe5zRJJ/Qc3uHPL99GA5vw9ColqPUPZKWustvXx39/kgmzn2fDe6/a9bkoN0l+OJHOpVwDI7qGoXPB1EJxfXBNLlfmFZfx3bFULuaW8OALv7um3vXYAof7WJYPAX757gvKjX7AY8B5IBzkGRBPc0OfCwSHa9izaQ0n9k62O9Jw9HR3NDIxt5d07lVC514l3DkLfttfwIGtARzZdRup527j09fzGTYph4j25qmNJc7if/+cR/r5M7XGZ5xIySejsJSx0REEqeAshRNcczaGoxfz+Dkho1rUY30oytfw3Sde7N4YgqlCD2xG6/EBMTfDkZ3f2rU5CI2GJ5dusI40qiqZmqx4cQ4BwaHVRia1tc9O8+Dn9UH8srkZxlINQvwfUj4NnLTbvjYjpaeHhtHdw+kY5mf3c8W1T6MaH4UQq4FbgOZAGrBASvlfR+2vhGIoNpaz9UQ6Z9LrTn5hLzzaUKRh29ogtn8ZiLFEQ3D4HrJS/oSH5xkqyowMvH0yo6bNtuvA9OhrK2jV4QaXnk9NWe9+/F0O7Yhi29pAjAaBEB8gTX8DMgGzQogZPKpOI6UQMKBDCAM6hFwReRXuTaO6REspp0gpI6SUOillZG1K4UpwJqOQ/8Wdc0opQPVwZ1MFxH0TwN+nt2PLpyHc0LeYef8+R8sOCxk0ridzF3/OTWOnUJCT6dCBac/GNXUe83KxyLrzq8WMvi+b51eeI7zND0jTw8AZYB4Ic7yFM0ZKKWHPmSw2HU6hTLlTKxzQpIOojOUmfk7I4OjFPKfa2/oIDAIWA33w8jnEH5+Hrv28He6fn5XO64+Oo6QwH2myvamc9TW4PFnNaHWedOs/FA9dDCd/vZPi/MF4eifQte9HmEz7ap2W1CQswItxsS0J0Ovqbqy4Jrjmg6guZBfzv7hzTisF+D3c2cOzJfAhsBMhwmlzwz8pLe7FsT2v1br/96veo6Qgjz7Dx9uETdfHdbo+slqOITTmZC29bhnD9AVLuP+vj/DKujAefOESeu8OHN71CudP3k/SiQSWPHUfF8+cqNNbMz2/lDV7z5OSV+KwjeL6pMmNGIzlJnaezuBwcp49T+NakRL+/cx3nIq/FwgC3gReAYqrtXPWM9LS1mKDcCZrU31Yt3gBuzfZn6ZUlbGkUMM7T5wn7fxwtB5JVJTfS4u22aSfP+OUXB4awege4XRu4e9S+RXuxzU5YriQXcwncec4dMFWKVQNabZHTroHH8xvxan4xzA7JfUB/kpNpWDvyW8vsUqz5uH0HTmhmg3C1RTkZtFv5ARu6DcEjVZrPXZVGf8yNobnJ3Yi7fwIYBgV5QLYSdq5SUgp6swmBWZ/h81HUvglMcvl56BomjQJP4ZiYznbEzI5kZLvsI2jfItSwv4tAWx4LxSTSXDnrHRiBhvZuLwtR3afcirU2Z7RsfuAYdbj3PXYArQaQai/F828dfh6eeDnpcXXywOtEEjMeRRMJigpKyevpMz8V1xGXkk5JgdDH4u9YN3iBUiTya6vxPyPtvL1stc4vGsL5cafgFjg38DfgVFodQ8Te3M3xj/yTLW+a67MSAm7z2SRU1zGyG4t0GpUXMT1jFsrBiklh5Pz2H0mC0OZOaS45g+6thwK8z8+xtrFLTi2x48OPYqZ8nQaIRFlQPVQZylNxAwehV+zYKc9IwtzMujUwo+Wgd5ENNMT6ud1WUFLpeUVXMo1kJxTTHJOCWn5BpvRkD2vTAsWpVVRZkRoNEhTHjAF+AZYQkXZXkoKllsVieX6Vc0NUVWRnkjJp8BQxrjYliqd/XWMW9sY1h1I5kJ29aH+uncWVot0dBQA1bHna2xa3pnSEsHtMzK5eUIumir3bX0djQC8dBo6NPejY5gfbUN8roiLcWFpOSdT8zmRUkBGgX27Rk0s52KJyyjIzcQ/sDmePrFkJr9GcUEHBow5Sdr5B0g6vs+pFZVgX0/u7NmKZj5qxeJaotGjK+uDs4ph2fYzFJWaRwq1hSP3HzmhSr5FHWGtvyb9wjAiOxmY9kwqLdqYvRWdrftQk7AAL2IjA+kS7t+o8QZZhaUcvpjHsYt5lFWYv6f6nsNn/3qVX765BZgOfA9MBarbEnoPH2fXMcrHU8u42Ja0DHS8hKtoWlxzxsfaMjFbhtr3/L/v0fueIv3CUEZMzuLxt89blQJUt0PUhRDQuYU/9/ZrzbQb29KjVbNGD0IK8fNiWJcwHhrcgQEdQvD21Dp9DpYU9b988zEwA3gIGAL8CpirLAmNBiGEQ8eoYmMFXxxIJiGtwNWnpnBz3HbEkJKSwtDb7mDac78/GR1lYp7454X88Fkw330cQrPm5Uz9SypRMb+vzdcn+YkQ0CnMnwEdgglpYC4DV+Js0hkLlimW2ShpOfdewJdAKDCd2CF5VrtKbdMoIWBQx+b0axfssI2iadDkRwwvv/wypw9XfzLaz8S8g3m3pfDNyubEDilg3r/PVVMKYH+50d6yZMcwP+4b0JbbYyLcSikAJCYmMnXqVHx8zElfvfTe9KnFqaqmUdLMQaAvvgFngDVcSryfCX9eYE1k62i5V0rYeSqTrcfTMNVRl1NxbeB2isHb2xshBO+//75N4tTpC5awYNXPVW7yicAhNNp+3Dk7gfueS8Xbz2xYq/pDrz3cGUL9vbi7TyTjYls2OOPRlSIiIoKAgAAMBgN6vZ4yYykxHSKI7tTW4T6WKVZUTH/C23aiY88BDBo3knbd/8KNt+WRkTyVFQtbYigWNlMUe4riyMU8voy/aF0hspdkV3Ft4HbLlYmJidXSsVWNYATzk1DnGUJZ6bvAw8BeYoesZsids6r1U9Ovwd6Sn7enlpuiQohu1axJ5DNMS0tj5syZPPLIIyxbtoyUlBQm9o7k+KV8fk7IsN6wFmqbHkiZRssOpWxYGsxf7ywGdgLSutxr8euouZx5LquYz/df4I7YVtWS7L733ntX5qQVVwW3tDHMmjWLZcuWofXQUVajqtOFBC+WztNgNLSi/+gzaD1epjA31XoTOGtP6N4ygCGdQ6+ZtfpiYzk/nczgZGr9DIW//mjk0zfaY6ooAMYD++y2q3r9HF1jvV5PSYmKu3BnmrSNIS0tDSGEVSkA7N64hidHLedfcyLw9gtn1usXmfyUZPR9j1CUn2sd8tZlTwj00XF3n0hGdQ+/ZpQCgI+nB2OiIxgXG1Gv8+o93JOYQQuBEuBn4C48PL1qra1R8xr7+Pgwbdo0zp4967oTUjiNs/4u9cEtFcP69eu5cOEC/W8dV/nji0SIH4DX6do/n3n/PkennuYnU825sSN7QmBIGP3aBXP/gLa0Dr52qzZ1DPPnvgFt6nWOFRWH6DfqdSI6lADrKDc+YdceY7E7ULnEaWlTUmLA39//sgvKKi4PKSX7krJZvfe8y2t/uJ2NwUJERAR6Hz/KSicAS5FSR1TMMh5++RaEsB3OWubGWp0nXt4+9Ln1ToZOfJC4zZ9RkpfFPf0iiWh2fTjq+Ot13NW7FQfO5bD7TBYVdawkTF+wpPJ6rgBWAP8Eoig3zkZohDVArKoSrmmziU9IorS8Ai+Pa2cU5s7kFhv5/lgaF3OvzNTNLW0MADk50KHLNnIzhhHRPpuI9q9QVnrUaktw5Aqt0Wo5sPUrBt4+mUlzFxLbOpDBHZtftxmSU/JK2HQ4hQJDea3tfvd7+IFy41+B+fgHHeTPb5Xy5sxbHCaNeWPTEev7ED9Pxse2JNBHJZy9Uljih3aeziQzLdXqBfvCPTc5ZUBvVBuDEOIPQoiTQojTQohnXdHn1KmQn30Ltz2QyZPvZXLfszOrWdlrThnKSg38+uP/VasP+f9GdmFMz7bXrVIAiGjmzbQb29KmjqnF734PBjw8XwFmUJgbzcqXejHnrR0Ok8ZUJavQyKd7z5OYUT3FnlrWdA05RUbWHkjmx9/SMZab6uXJW18aPGIQQmiBBGAkkIzZrD1FSnnc0T7OjBiOHIHVu5MJ7VDssE3VQKjt61eScHAPxQW5lJUa8Pb2ZuLEibz55ptq7ov5SbPnTBZ7k7IdJripGVh2KbEVKUn/QOclad/tBQ7v/Ifd/Wqu+AgB/dsFMzDKXNZu9uzZ/Oc//+HRRx9Vy5qXgclkrrcal5hFWYVs0KpQowVRCSEGAgullKMr3z8HIKW0/yvi8oKonOGLdxaye9NneHl5YjQa1Q/RDmcyCvn2aKrT6fVTz3nywfxW5GZUEBX9Lzw8vyXh113Wwru1Fdh5ZmxMtZUlC2pZ03lS8wz88Fsa6fm/X0d70+hJd03krbfeqvMh2JhTiVbAhSrvkyu31RToESHEfiHE/oyMuqtG15cQP08CRTGzZs0kLi6OmTNnqqGrHaJC/ZjcrzWBToZTh7c18sQ754nsBGcOP4uh6CFMFfaTxtTk+Y+20v/W8ei9zUZftazpPKXlFWz7LZ01+85XUwpgf+UtICDApSNjV6xK2LN42AxDrmSJum4tAxh+Qxh//PpL67alS10/77pWCPHzYkr/Nmw6nML5bMdTNQv+QRX8+Y0LrHotnCO7/kh4uwimzDOw97s1tSabDQgJw8PLh1KDAU8vLwwGg8t/wNciv6XmsyMhk8Ji3Z1JAAAV6ElEQVRSxwbjqqtC29evZO3atSxYsMBl19YVI4ZkoHWV95HAJRf0Wyc6rWBU9xaM7h5+XRsYLwe9TsuEXq3o2ca5greeeskDf0th2KRsUpNGsnnlcMbMeNFuAFbV95Yf8GNvf86tE6eRfKlRfhpNkvQCA5/vv8A3R1JrVQr5WekU5ecyctpsWkXdgM5LT05ODi+99JLLZHGFjcEDs/FxBHARs/FxqpTymKN9LsfGUDNBSYifJ7dHu18UZFMk/kIuP5/McJh7siZx3wSw7p0WhEYaefili2xbO79aVq2aWbaqotdpuaVLKF0jAq7AmTRNSowV7D6TydGL+U59B5brixB2s3HVZsNp7BJ1Y4C3AS3woZTy1draX45iqPpje/7vb3Fr1xZ4eqhRQkNJSUlh8uTJvPneCvamm5w2Sp466M37z/gBZcAEYLfDtvZyRnQI9WVE1xb4ebmtj90Vp7zCRPyFXPYmZVNa9vt1d5Sly3EZg2ZAGD4+F5kwYUKtK3GNXaJus5Sys5Qyqi6lUF8smYh2b1xt9U+4PaYlzfx9XXmY6xZLhOSKpW9yb7/WBHhXN0o6ytPQqVcJc95KxMu7HNgGTMfD04ugFq1qjbOwkJhRxCd7znH8kuPM39cqUkqOX8pn5e4kdpzKrKYUwHGmMdssZt74Bc4DEhBiNSUlrrPhuP0j13IxPFXAjkupmvfCZDLx/vvvE+qvZ87IbkQ001vb1eZE0yHan9ibXwB+Aj6k3Pg6Ok9zchiL09np+DiHMhjKKvjuWCpr918gs9D1gUDuyOn0Av73y3m+O5Zq441q7yH45KguPH17tE2MitajP+XGrRTmvoF/UC6rV/sxa5brVuLcXjEEhIQREhxIeZkRvV6vLNsuomZGqKoK964+kTzj4Edas3BNSdF5brr9E/qMOAM8TlbKf+g36mHmLv6c8LYdyc/OqNMzLzmnhFVx5/npZLpNTolrhTMZhaz65Rz/dyiFTAfRkI4ig3vdMoazR/ezcfmbHNy+h+aRGzFVxKHz7E5kxzdYsBruuaczS5cuZf369S6R1+0neH3aBrG9osgmQYmiYdTMCFVT4SadTeSBRx/j5y2bq8Wi1Cxc87ubegWdeqWydvHN7P2uLXu/mwicBqrX+nBU9NckJQfP55KQVkD/9ubkOU296I3JJDmdUcj+pBzS8m3zddbEkZu/GcH+rc2AvZQUhDBofB5jHszE2++OKyK7WyuG0d3DaRviy5ANv2tB5Z/gOuxlhLLQsmVLOkaG8YPR6JQzE0D/Ufm0aGNkxYvhFOTEodE8TkX5fxwqFXsUlZode/YnZTOgQwjdIgLQNDEFUVZh4tilfH49l0NeSVm99q3pn5BwcA95mS2ApZirs+8GRrHr60P88q3rqqvXxG2jKxVXn4kTJxIREcHtk+7nn28vITczvc6iPACFeRremlVCXmY0QrMCaZrNTWMnXFbR30AfHX3aBtE1IsDtfVWyi4wcTs7lREqBS6ZERfkalj6VSOq5EZhrgfwF+Bidl5eNK/oTt3ZyaXSlW48YFFeXqvPVfn168VX8xTrDtwH8mplo3flp/INmknxqOj4Bw8m4+PJlyZBbXMYPJ9LZdTqL7i0DiG0dSDNv96mOZSw3kZhZyJHkPJJzXBP/YaqAX75txuYVzSnKb09Eu60EhS/neNw6hNA4NXprKO6tghVuQ6i/F5P7tyEswDmHshkL3+HJpd14+JWLaEQrko4vY9vngnefvI+LZ07UWpnc3hKpoayCA+dyWLkriQ0Hkzl6Me+qGSrLK0zsPnKK6L4DefPLOL45ksrx0+dqPSdnSTyq5+3H27B2cQtatCnlqfcv8PSy9mi1eQwaN5Un39twxaqrV0VNJRT1wlhu4pujKSRmFDm9T36Wlk9fDyfhoC+wltDIf5J58SB9RtxBdtpFG0ee2jwnq+KhEbQJ8SEq1I/WQT5XtM5mVmEp57KLuZBtLj786aIXnPL2dLakYHaaBxuW+nMsLhT/4FLGP5JN72EFOJu83NVTCaUYFPVGSsmOU5kcOJdTbXvtHntGYB7wKpAOPAJsBsyl8t769kS9KobZw1/vQWSQDxHN9AT7ehLs64lvPT0rpZQUlpaTXlBKen4p6QUG0vNLrbELjr0P7ctcl5IzFAu2fR7MtnVBmMrLMJn+yY1/SOLeJ5+rl9xKMSjchiPJefz4W7rVv7+2p+bv5fK6ASuBGOAj4AkgFwCNh46eN4+2SddXM99DfQr76nVamnnr0Os0eHlo0es06LQaKqTEZJJUmCTlJkmhoZyC0nKKSstrzZFZMxeCh6cX/kHNKcjJpNxYapX58I7vKC8z2uxvURgV5bBtrZZvPw7CVBEMfAo8gzkmsXpbZ1DGR4XbEB3ZjEAfHR0igu0m5rX8sKuXyzuENPUF5gN/BUaCeBzkF/QedjuennqHFcMs1CwmVBuGsgqX2iJq+hpUlBnx9PK2entaZJ7/8Q92c5KOffgZDm33Y/OK5mRc9AS20bXfZrz9fuPI7kzKSqnX8u6VQikGRYNoHezDkeMnmTHrcfZv/96hM5RlfT7tQiKFOVkUFy4lP+tL4EOQ64Dv2b/lceAkQqNh7uLP2b5+JfHbv2HUfX8mIDjUYWbw+jxZXUHNDNlHdm2xqXJm66xUSnH+EP77tz5cPKMHjmJefvyGE1Vq/DjrM1IVH0+tyyupKcWgaDA3RLWlR/sI4rY6ftLX9H9Y8eIc9D5+5Oc8ScKBaKR8EThCWJv/46EXWxLaKhidl56Sgjw2Ln+T7LSLzF38Gds+X15tGF+1fGFjUfVc7npsAXc9tqDaezBPOeJ3fEvvEXfSpvNf+H5VGL/t70DzlkYmzD5J0vHnObpnG2WlZmUgTSZ6DLqVWyfPtCqXuvD21NKnbRCxkc7l1KgPSjEoXEJmRjozZ86k/233sPyDD+r8YVturnWLFwDvotWtp6LsRdLPz+Af0zMw17bYAEj2bzVn5lr054kMvO2ealONQmMpezaucdp5ypF9oj52C2f6DwxrQ3H+WE788ioHtrYhONzIuD+l0mdEPlqtIO2cT7XzAPD1D6RV1A3VFI09vHQaercJolebwCtWx0P5MShcwvr163nvvfd4cNwtrFz+H2a+4lwSXsuw/Il3ljJo3A9ExcwhKsYDeAs4A8wGzHUqpMnE7k1rkFLaTCnsBXjZw1G0qKNq3xafC0e+FzV9LhZOHU/ikUH8+sMi4BOKC/KBB8jLCqb/qHy02t/PGyFqPY+afWs1gp5tApl+U3sGdAi5osV91KqE4oqQU2Rk4+FLZBbaWuad4b8vrOZY3DBgCJACvI+H50fEDI5l2KSHqk0p6spWDc4vM1qwVPsOaxNF+vkz1v81V1zWvbOQ3RvXgLgR5J+AyYAP5lD014Fv6D18nF3ZHBVNsrStusrz11ffYlDHkAYX81GrEoqrSpCvJ5P7t2Hr8TR+q2cFbgCNdhc3jU0gJelbzh69GXiJcuPzXEiIo8zYHi9v2/qklvqa9qYE8z/aavcmrKlkLFhSpqWdO13tv8XgaaY1cA9wAGQvoBD4GPgPEF+pXIRDQ6KjOquv/HGEzUhi7MbVjZp2X00lFFcMnVbDbdERjOzWAp22flbz6QuWcPfjC/Brtp9B4/7Hgwt+JLzdz2Rd6s87c9uw/4dFtIz6kilPf8vA26fYra9ZFUc3YauortW2g9kYqNXZfzJrPLrhqX+WFm3SgPPAm0A5Go/HCA4fgBCzEcK8QhIzeJTVfdlRJizLVGru4s+tbed/tJUbr3La/QZNJYQQk4CFQFegv5TSqfmBmkpcf2QXGdl8JKXBJdtLijQc2eXHwZ/8OfWrDyaTICTCSHba/5CmHzCnmfs9fLzqUmbNymXHftnGrNdW8v4zD9LtxmEMnfggq/75FKmVowOEBmQboD/mQmsjgbYA+AQk0TwijvMn5+PhmUxFmZGgFq3o2m9ItWVLq5G1yrRg1LTZDg2dnh4a+rcPZvk//8ryDz7A09O1xZMaxfNRCNEVMGEeO81TikFRG+UVJnaeziT+Qq7DMnn1oTBXy+Gdfvy2z5fTh/UYiiwz47MIzW+EtjIw4LYeRHbS4xNQgY+fCR//CnReknXvvMieTRsIbd2FjAtpICJBtsRc/aAD0AvoCQRV9pkL/AhsAb4HEgGzLeLJpRtsFIGF2mwbN42dYrVXaISgR6sABkaF4OPpYQ15r5orwxXZmRo7S/RPKMWgcJIL2cVsOZ5W7yQmtS0pmirg41dXcHhnGULTH2nqjNB0Q5oux1hXghDHCInIZMCYjnSM1eEfdIFN/3VsKKxN5qq2DXvoPL24lJVH80YoheB2xkchxCOYI2do06ZNYx1W4Ya0DvbhvgFt2XU6k0PJv48e6vIlqM0VWqMFKQ8waFwoA8boidv8BnlZWYx/5H2y0zwoKdBSXKAlJ72Yo7t/IT35JKaKEsAIlKDxSCcg0EBu5gG0ujxM5aV07j2Z4fcsBEoB+zaKunweqto2tDpPKsqMaDRaTKYKPL30jBl3B++/+3ajKIX6UKdiEEJsBexlXn1eSvmVswe6kiXqFE0PTw8Nw24Io1MLP378LZ2sQqPDG99ZV+iaHolmymjesqya0inO/4q0858hNOaCLUKjQVZINB6tGDRuRDUbQVVqukI7m3uh6n4WG4aXl56yMiMRocFumdhYTSUUVx1vb28MBtthtuXGr2u93x4WRTBh9vNseO9VgsNbcWDrVwy8fTIFOZkEBIda4zb8gkJo0bqDXRuBKwn00bHq1cfp3K41jz76qEttB87idlMJhcIRiYmJzJs3jw0bNlBSUmIThOVoqbG2Ybxl9PHW7AkgJYlHzJFKFh+Exgy8auato3/7YLpFBDB98/9Zt7tzYuMGKQYhxATgXSAU2CSEiJdSjnaJZIrrBksq+9LSUvR6PUajkbCQwGo3vrPDeGc8HC2eiFeasAAv+rQNonOYf5PLdK1cohVugb3luaUrVrH7dBbns4ud7qd6UpjqCsLi5lxXyriGIAS0DfGhd5sg2oa4XxlFNZVQNCmqzrOrDrHv6hNJWr7BWoymtuxKQI2kMJpq1aBjBo/Cr1lwgxO22sPHU0v3ls3o0SqgwfEM7oBSDAq3p0WAnj/0COfmTs05cjGPoxfzak1jXzMpTFXjYl0hzfXBQyNo19yXLuH+RIX6NfnKWVVRUwlFk0NKSUqegZNpBZxKK6CotPHSyHt6aGgd7EOnMD86hPpe0dDnK4GaSiiuWYQQtAz0pmWgN7d0DuVQwllm/PF+/vzyuxh0ARjLTXV34iSeHhpC/b1oHeRDmxAfIgL0Tc6QeDkoxaBo0gghWLb4TQ7tj+PAl8t5d8lSsgpLySoyklNkJKvISF5JGYayCkrLTXaVhodG4Kf3wNfLA38vDwK8dTT38yLM34tAH53L8yk2BdRUQtFkceQYVVveggqTpKzChEYItBqBRnBd3fjOTiVUPgZFkyUxMZGpU6fi4+MDOJe3QKsR6HVaPD00aDXiulIK9UEpBkWTxeIYZTAY0Ov1GAwGAgIC3DL2oKmhFIOiSZOWlsbMmTOJi4tj5syZpKamXm2RrgmUjUGhuI5QNgaFQnHZKMWgUChsUIpBoVDYoBSDQqGwQSkGhUJhg1IMCoXCBqUYFAqFDUoxKBQKGxqkGIQQbwghfhNCHBZCbBBCBLpKMIVCcfVo6IhhC9BDShkDJADPNVwkhUJxtWmQYpBSfi+ltOTYisNc+E+hUDRxXGljmAF848L+FArFVcIlJeqEEM8D5cCqWvpRtSsViiZCnYpBSnlrbZ8LIR4AxgIjZC2hmqp2pULRdGhoJao/AM8AQ6WUzlcFUSgUbk1DbQxLAH9gixAiXgjxbxfIpFAorjINGjFIKTu6ShCFQuE+KM9HhUJhg1IMCoXCBqUYFAqFDUoxKBQKG5RiUCgUNijFoFAobFCKQaFQ2KAUg0KhsEEpBoVCYYNSDAqFwgalGBQKhQ1KMSgUChuUYlAoFDYoxaBQKGxQikGhUNigFINCobBBKQaFQmGDUgwKhcKGhpaoe7myPF28EOJ7IURLVwmmUCiuHg0dMbwhpYyRUvYENgIvuEAmhUJxlWloibr8Km99AVUvQqG4BmhQlmgAIcSrwB+BPGBYgyVSKBRXHVFL8ShzAydK1FW2ew7QSykXOOjHWqIO6AKcdEK+5kCmE+2uJu4uo7vLB+4vo7vLB87L2FZKGVpXozoVg7MIIdoCm6SUPVzSobnP/VLKvq7q70rg7jK6u3zg/jK6u3zgehkbuirRqcrb8cBvDRNHoVC4Aw21MfxTCNEFMAHngJkNF0mhUFxtGlqi7i5XCeKAZVe4f1fg7jK6u3zg/jK6u3zgYhldZmNQKBTXDsolWqFQ2OAWikEI8QchxEkhxGkhxLN2PvcSQnxW+fkvQoh2bibfk0KI45Xu4T9UrtA0KnXJWKXd3UIIKYRodCu7MzIKIe6pvJbHhBCfupN8Qog2QohtQoiDld/1mEaW70MhRLoQ4qiDz4UQ4p1K+Q8LIXpf9sGklFf1D9ACZ4AOgCdwCOhWo81s4N+VrycDn7mZfMMAn8rXsxpTPmdlrGznD2wH4oC+7iYj0Ak4CARVvg9zM/mWAbMqX3cDkhr5Gg4BegNHHXw+BvgGEMAA4JfLPZY7jBj6A6ellIlSSiOwBrijRps7gI8qX68DRgghhLvIJ6XcJqUsrnwbB0Q2kmxOy1jJy8DrgKExhavEGRn/BCyVUuYASCnT3Uw+CQRUvm4GXGpE+ZBSbgeya2lyB/CxNBMHBAohIi7nWO6gGFoBF6q8T67cZreNlLIcs/t1SKNI55x8VXkIs9ZuTOqUUQjRC2gtpdzYmIJVwZnr2BnoLITYJYSIE0L8odGkc06+hcB9QohkYDPwWOOI5jT1/a06pMGxEi7A3pO/5lKJM22uFE4fWwhxH9AXGHpFJbJzaDvbrDIKITTAv4AHG0sgOzhzHT0wTyduwTzq2iGE6CGlzL3CsoFz8k0BVkop3xJCDAQ+qZTPdOXFcwqX3SfuMGJIBlpXeR+J7RDN2kYI4YF5GFfbkMqVOCMfQohbgeeB8VLK0kaSzUJdMvoDPYCfhBBJmOefXzeyAdLZ7/krKWWZlPIs5niaTjQOzsj3EPA5gJRyD6DHHKPgLjj1W3WKxjSeODCYeACJQHt+N/p0r9Hmz1Q3Pn7uZvL1wmy46uSu17BG+59ofOOjM9fxD8BHla+bYx4Wh7iRfN8AD1a+7lp504lGvo7tcGx8vJ3qxse9l32cxjypWk52DJBQeXM9X7ntJcxPXzBr5rXAaWAv0MHN5NsKpAHxlX9fu9s1rNG20RWDk9dRAIuA48ARYLKbydcN2FWpNOKBUY0s32ogBSjDPDp4CHMYwswq129ppfxHGvIdK89HhUJhgzvYGBQKhZuhFINCobBBKQaFQmGDUgwKhcIGpRgUCoUNSjEoFAoblGJQKBQ2KMWgUChs+P+fWr/AqVxL/AAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 288x216 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "with torch.no_grad():\n",
+    "    # Initialize plot\n",
+    "    f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
+    "\n",
+    "    # Get upper and lower confidence bounds\n",
+    "    lower, upper = observed_pred.confidence_region()\n",
+    "    # Plot training data as black stars\n",
+    "    ax.plot(train_x.cpu().numpy(), train_y.cpu().numpy(), 'k*')\n",
+    "    # Plot predictive means as blue line\n",
+    "    ax.plot(test_x.cpu().numpy(), observed_pred.mean.cpu().numpy(), 'b')\n",
+    "    # Shade between the lower and upper confidence bounds\n",
+    "    ax.fill_between(test_x.squeeze().cpu().numpy(), lower.cpu().numpy(), upper.cpu().numpy(), alpha=0.5)\n",
+    "    ax.set_ylim([-3, 3])\n",
+    "    ax.legend(['Observed Data', 'Mean', 'Confidence'])\n",
+    "    plt.savefig('res.pdf')"
+   ]
+  }
+ ],
+ "metadata": {
+  "anaconda-cloud": {},
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}
diff --git a/examples/01_Simple_GP_Regression/index.rst b/examples/01_Simple_GP_Regression/index.rst
index 1bb65c7a..a4e94b45 100644
--- a/examples/01_Simple_GP_Regression/index.rst
+++ b/examples/01_Simple_GP_Regression/index.rst
@@ -5,6 +5,7 @@
    :hidden:
 
    Simple_GP_Regression.ipynb
+   Simple_MultiGPU_GP_Regression.ipynb
    Simple_GP_Regression_With_LOVE_Fast_Variances_CUDA.ipynb
    Simple_Batch_Mode_GP_Regression.ipynb
    Spectral_Mixture_GP_Regression.ipynb
diff --git a/gpytorch/kernels/__init__.py b/gpytorch/kernels/__init__.py
index 591786d3..1610a984 100644
--- a/gpytorch/kernels/__init__.py
+++ b/gpytorch/kernels/__init__.py
@@ -2,6 +2,7 @@
 
 from .additive_structure_kernel import AdditiveStructureKernel
 from .cosine_kernel import CosineKernel
+from .multi_device_kernel import MultiDeviceKernel
 from .grid_interpolation_kernel import GridInterpolationKernel
 from .grid_kernel import GridKernel
 from .index_kernel import IndexKernel
@@ -23,6 +24,7 @@ __all__ = [
     "Kernel",
     "AdditiveKernel",
     "AdditiveStructureKernel",
+    "MultiDeviceKernel",
     "CosineKernel",
     "GridKernel",
     "GridInterpolationKernel",
diff --git a/gpytorch/kernels/multi_device_kernel.py b/gpytorch/kernels/multi_device_kernel.py
new file mode 100644
index 00000000..96b40405
--- /dev/null
+++ b/gpytorch/kernels/multi_device_kernel.py
@@ -0,0 +1,62 @@
+#!/usr/bin/env python3
+
+import torch
+from torch.nn.parallel import DataParallel
+from .kernel import Kernel
+from ..lazy import CatLazyTensor, lazify
+from .. import settings
+
+
+class MultiDeviceKernel(DataParallel, Kernel):
+    r"""
+    Allocates the covariance matrix on distributed devices, e.g. multiple GPUs.
+
+    Args:
+        - :attr:`base_kernel`: Base kernel to distribute
+        - :attr:`device_ids`: list of `torch.device` objects to place kernel chunks on
+        - :attr:`output_device`: Device where outputs will be placed
+    """
+
+    def __init__(self, base_kernel, device_ids, output_device=None, **kwargs):
+        DataParallel.__init__(self,
+                              module=base_kernel,
+                              device_ids=device_ids,
+                              output_device=output_device,
+                              dim=-2)
+
+        self.output_device = output_device if output_device else device_ids[0]
+
+        self.__cached_x1 = torch.empty(1)
+        self.__cached_x2 = torch.empty(1)
+
+    def forward(self, x1, x2, diag=False, **kwargs):
+        if diag:
+            return self.module.forward(x1, x2, diag=True, **kwargs).to(self.output_device)
+
+        if not x1.device == self.__cached_x1.device or not torch.equal(x1, self.__cached_x1):
+            self._x1_scattered, self._kwargs = self.scatter((x1,), kwargs, self.device_ids)
+            self.__cached_x1 = x1
+
+        if not x2.device == self.__cached_x2.device or not torch.equal(x2, self.__cached_x2):
+            self._x2_subs = [x2.to(x1_[0].device) for x1_ in self._x1_scattered]
+            self.__cached_x2 = x2
+
+        inputs = tuple((x1_[0], x2_) for x1_, x2_ in zip(self._x1_scattered, self._x2_subs))
+
+        if not self.device_ids:
+            return self.module.forward(*inputs, **self._kwargs)
+
+        if len(self.device_ids) == 1:
+            return self.module.forward(*inputs[0], **self._kwargs[0])
+
+        # Can't cache the replication because the base kernel module can change every time (e.g. param updates)
+        replicas = self.replicate(self.module, self.device_ids[:len(inputs)])
+
+        # TODO: parallel_apply might be too heavyweight in some cases?
+        with settings.lazily_evaluate_kernels(False):
+            outputs = self.parallel_apply(replicas, inputs, self._kwargs)
+
+        return self.gather(outputs, self.output_device)
+
+    def gather(self, outputs, output_device):
+        return CatLazyTensor(*[lazify(o) for o in outputs], dim=self.dim, output_device=self.output_device)
diff --git a/gpytorch/lazy/__init__.py b/gpytorch/lazy/__init__.py
index e4bddb07..b74bda3e 100644
--- a/gpytorch/lazy/__init__.py
+++ b/gpytorch/lazy/__init__.py
@@ -6,6 +6,7 @@ from .batch_repeat_lazy_tensor import BatchRepeatLazyTensor
 from .block_lazy_tensor import BlockLazyTensor
 from .block_diag_lazy_tensor import BlockDiagLazyTensor
 from .cached_cg_lazy_tensor import CachedCGLazyTensor
+from .cat_lazy_tensor import CatLazyTensor
 from .chol_lazy_tensor import CholLazyTensor
 from .constant_mul_lazy_tensor import ConstantMulLazyTensor
 from .diag_lazy_tensor import DiagLazyTensor
@@ -33,6 +34,7 @@ __all__ = [
     "BlockLazyTensor",
     "BlockDiagLazyTensor",
     "CachedCGLazyTensor",
+    "CatLazyTensor",
     "CholLazyTensor",
     "ConstantMulLazyTensor",
     "DiagLazyTensor",
diff --git a/gpytorch/lazy/cat_lazy_tensor.py b/gpytorch/lazy/cat_lazy_tensor.py
new file mode 100644
index 00000000..968f361a
--- /dev/null
+++ b/gpytorch/lazy/cat_lazy_tensor.py
@@ -0,0 +1,346 @@
+#!/usr/bin/env python3
+
+import torch
+import itertools
+from collections import defaultdict
+from .lazy_tensor import LazyTensor, delazify
+from .non_lazy_tensor import lazify
+
+
+class CatLazyTensor(LazyTensor):
+    r"""
+    A `LazyTensor` that represents the concatenation of other lazy tensors.
+    Each LazyTensor must have the same shape except in the concatenating
+    dimension.
+
+    Args:
+        - :attr:`lazy_tensors` (list of LazyTensors):
+            A list of LazyTensors whose sizes are the same except in
+            concatenating dimension :attr:`dim`
+        - :attr:`dim` (int):
+            The concatenating dimension which can be a batch dimension.
+        - :attr:`output_device` (torch.device):
+            The CatLazyTensor will appear to appear on :attr:`output_device`
+            and place any output `torch.Tensors` on :attr:`output_device`
+    """
+
+    def __init__(self, *lazy_tensors, dim=0, output_device=None):
+        if len(lazy_tensors) == 0:
+            raise RuntimeError("List of LazyTensors must be non-empty")
+        elif len(lazy_tensors) == 1:
+            raise RuntimeError("Why are we trying to concatenate a single LazyTensor?")
+        if not all([isinstance(t, LazyTensor) for t in lazy_tensors]):
+            raise RuntimeError("CatLazyTensor requires a list of all LazyTensors")
+
+        super().__init__(*lazy_tensors, dim=dim, output_device=output_device)
+
+        def remove_dim(tuple, dim):
+            return tuple[:dim] + tuple[dim + 1:]
+
+        rep_tensor = lazy_tensors[0]
+        ndims = rep_tensor.ndimension()
+        if dim < 0:
+            dim = ndims + dim
+        pre_cat_size = tuple(rep_tensor.size()[:dim])
+        post_cat_size = tuple(rep_tensor.size()[dim + 1:])
+
+        cat_dim_len = 0
+        cat_dim_sizes = []
+        tensor_idx_to_start_idx = []
+        for t in lazy_tensors:
+            if t.ndimension() != ndims:
+                raise RuntimeError("All tensors must have the same number of dimensions")
+            if remove_dim(t.size(), dim) != remove_dim(rep_tensor.size(), dim):
+                raise RuntimeError("All LazyTensors must have the same size in "
+                                   "the non-concatenation dimension")
+            tensor_idx_to_start_idx.append(cat_dim_len)
+            cat_dim_size = t.size()[dim]
+            cat_dim_len += cat_dim_size
+            cat_dim_sizes.append(cat_dim_size)
+
+        # using itertools to more quickly join list of lists
+        idx_to_tensor_idx = [[t_idx] * size for t_idx, size in enumerate(cat_dim_sizes)]
+        idx_to_tensor_idx = list(itertools.chain.from_iterable(idx_to_tensor_idx))
+
+        self.lazy_tensors = lazy_tensors
+        self.pre_cat_size = pre_cat_size
+        self.post_cat_size = post_cat_size
+        self.cat_dim_sizes = cat_dim_sizes
+        self.cat_dim_len = cat_dim_len
+        # can't call this attribute self.dim because LazyTensor has a dim() function
+        self.cat_dim = dim
+        self.idx_to_tensor_idx = idx_to_tensor_idx
+        self.tensor_idx_to_start_idx = tensor_idx_to_start_idx
+        self.tensor_idx_to_end_idx = tensor_idx_to_start_idx[1:] + [sum(self.cat_dim_sizes)]
+        self.output_device = output_device
+
+    def _split_slice(self, slice_idx):
+        """
+        Splits a slice(a, b, None) in to a list of slices [slice(a1, b1, None), slice(a2, b2, None), ...]
+        so that each slice in the list slices in to a single tensor that we have concatenated with this LazyTensor.
+        """
+        if slice_idx.step is not None:
+            # TODO: Add support for this eventually.
+            raise RuntimeError('Slicing a CatLazyTensor with a step is not currently supported!')
+        start_idx = slice_idx.start if slice_idx.start is not None else 0
+        stop_idx = slice_idx.stop if slice_idx.stop is not None else self.shape[self.cat_dim]
+
+        start_tensor_idx = self.idx_to_tensor_idx[start_idx]
+        stop_tensor_idx = self.idx_to_tensor_idx[stop_idx - 1]
+
+        if start_tensor_idx != stop_tensor_idx:
+            # By definition, stop is on a later tensor than start since they are in order.
+            end_idx = self.tensor_idx_to_end_idx[start_tensor_idx]
+            my_slice = slice(start_idx, end_idx)
+
+            if end_idx == stop_idx:
+                return [my_slice]
+            else:
+                # Keep splitting
+                return [my_slice] + self._split_slice(slice(end_idx, stop_idx, None))
+        else:
+            return [slice(start_idx, stop_idx, None)]
+
+    def _getitem(self, *indices):
+        squeeze = [isinstance(i, int) for i in indices]
+        indices = list(indices)
+        target_indices = indices[self.cat_dim]
+        new_cat_dim = (None if squeeze[self.cat_dim]
+                       else self.cat_dim - sum(squeeze[:self.cat_dim + 1]))
+
+        eval_result = (torch.is_tensor(indices[-2]) and torch.is_tensor(indices[-1]))
+        eval_result = eval_result or (isinstance(indices[-2], int) or isinstance(indices[-1], int))
+        maybe_lazify = delazify if eval_result else lazify
+
+        if new_cat_dim is None:
+            # target_indices must be a int so we can let the LazyTensor squeeze out cat_dim
+            t_idx = self.idx_to_tensor_idx[target_indices]
+            return maybe_lazify(self.lazy_tensors[t_idx]._getitem(*indices))
+
+        if all(torch.is_tensor(x) for x in indices):
+            left_indices, right_indices = indices[-2], indices[-1]
+            batch_indices = indices[:-2]
+            return maybe_lazify(self._get_indices(left_indices, right_indices, *batch_indices))
+
+        if isinstance(target_indices, slice):
+            if target_indices == slice(None, None, None):
+                res_list = [lazify(t._getitem(*indices)) for t in self.lazy_tensors]
+                return maybe_lazify(self.__class__(*res_list, dim=new_cat_dim, output_device=self.output_device))
+            else:
+                target_slices = self._split_slice(target_indices)
+                target_tensors = [self.idx_to_tensor_idx[sl.start] for sl in target_slices]
+
+                res_list = []
+                for idx, t_idx in zip(target_slices, target_tensors):
+                    shifted_start = idx.start - self.tensor_idx_to_start_idx[t_idx]
+                    shifted_stop = idx.stop - self.tensor_idx_to_start_idx[t_idx]
+                    shifted_slice = slice(shifted_start, shifted_stop, idx.step)
+                    indices[self.cat_dim] = shifted_slice
+                    res = lazify(self.lazy_tensors[t_idx]._getitem(*indices))
+                    res_list.append(res)
+                if len(res_list) == 1:
+                    result = res_list[0]
+                elif all([rl.dim() == 1 for rl in res_list]):
+                    return maybe_lazify(torch.cat([rl.evaluate().to(self.device) for rl in res_list]))
+                else:
+                    shape_diffs = torch.tensor(res_list[0].shape) - torch.tensor(res_list[1].shape)
+                    new_cat_dims = (shape_diffs != 0).nonzero()
+                    new_cat_dim = new_cat_dims.item() if new_cat_dims.numel() > 0 else self.cat_dim
+                    result = self.__class__(*res_list, dim=new_cat_dim, output_device=self.output_device)
+
+                return maybe_lazify(result.to(self.output_device))
+        elif torch.is_tensor(target_indices):
+            # this means another `indices` is a slice object
+            target_indices = [idx.item() for idx in target_indices]
+            target_tensors = [self.idx_to_tensor_idx[idx] for idx in target_indices]
+
+            res_list = []
+            curr_tensor, slice_indices = target_tensors[0], []
+            for idx, t_idx in zip(target_indices, target_tensors):
+                if t_idx != curr_tensor:
+                    indices[self.cat_dim] = torch.tensor(slice_indices)
+                    new_inds = [ind[:len(slice_indices)] if torch.is_tensor(ind) else ind for ind in indices]
+                    res = lazify(self.lazy_tensors[curr_tensor]._getitem(*new_inds))
+                    res_list.append(res)
+                    curr_tensor, slice_indices = t_idx, []
+                slice_indices.append(idx - self.tensor_idx_to_start_idx[t_idx])
+            indices[self.cat_dim] = torch.tensor(slice_indices)
+            new_inds = [ind[:len(slice_indices)] if torch.is_tensor(ind) else ind for ind in indices]
+            res = lazify(self.lazy_tensors[t_idx]._getitem(*new_inds))
+            res_list.append(res)
+
+            if len(res_list) == 1:
+                result = res_list[0]
+            else:
+                shape_diffs = torch.tensor(res_list[0].shape) - torch.tensor(res_list[1].shape)
+                new_cat_dims = (shape_diffs != 0).nonzero()
+                new_cat_dim = new_cat_dims.item() if new_cat_dims.numel() > 0 else self.cat_dim
+
+                result = self.__class__(*res_list, dim=new_cat_dim, output_device=self.output_device)
+
+            return maybe_lazify(result.to(self.output_device))
+
+    def _get_indices(self, left_indices, right_indices, *batch_indices):
+        # tensor indices must all have the same length
+        indices = list(batch_indices) + [left_indices, right_indices]
+        indices = torch.stack(indices, dim=0)
+        target_indices = indices[self.cat_dim, :]
+        target_tensors = [self.idx_to_tensor_idx[idx.item()] for idx in target_indices]
+        starting_indices = [self.tensor_idx_to_start_idx[t_idx] for t_idx in target_tensors]
+        local_indices = target_indices - torch.tensor(starting_indices)
+        indices[self.cat_dim, :] = local_indices
+        if len(set(target_tensors)) == 1:
+            # shortcut if target_indices are all on the same LazyTensor
+            left_indices, right_indices = indices[-2, :], indices[-1, :]
+            batch_indices = tuple(indices[:-2, :])
+            return self.lazy_tensors[target_tensors[0]]._get_indices(left_indices, right_indices, *batch_indices)
+
+        d = defaultdict(list)
+        for i, t_idx in enumerate(target_tensors):
+            d[t_idx].append(i)
+
+        res_list = []
+        for t_idx, slices in sorted(d.items()):
+            indices_ = indices[:, slices]
+            left_indices, right_indices = indices_[-2, :], indices_[-1, :]
+            batch_indices = tuple(indices_[:-2, :])
+            res = self.lazy_tensors[t_idx]._get_indices(left_indices,
+                                                        right_indices,
+                                                        *batch_indices)
+            res_list.append(res)
+        # collect all the res in res_list onto one device
+        res = torch.cat([r.to(self.device) for r in res_list], dim=0)
+
+        t_idx_to_res_idx = []
+        curr_idx = 0
+        for t_idx in sorted(d.keys()):
+            t_idx_to_res_idx.append(curr_idx)
+            curr_idx += len(d[t_idx])
+        lookup = []
+        # use the fact that order of elements retrieved from each LazyTensor is
+        # the same as the order they appear in target_indices
+        for t_idx in target_tensors:
+            idx = t_idx_to_res_idx[t_idx]
+            lookup.append(idx)
+            t_idx_to_res_idx[t_idx] += 1
+        return res[lookup]
+
+    def _matmul(self, rhs):
+        isvector = rhs.ndimension() == 1
+        if isvector:
+            rhs = rhs.unsqueeze(1)
+
+        output_device = (self.device if self.device is not None
+                         else rhs.device)
+        # make a copy of `rhs` on each device
+        rhs_ = []
+        for d in self.devices:
+            if d != rhs.device:
+                rhs_.append(rhs.to(d))
+            else:
+                rhs_.append(rhs)
+
+        if self.cat_dim == self.ndimension() - 2:
+            res_list = [t._matmul(rhs)
+                        for t, rhs in zip(self.lazy_tensors, rhs_)]
+            # copy result back to output device
+            res_list = [x.to(output_device) for x in res_list]
+            res = torch.cat(res_list, dim=self.cat_dim)
+        elif self.cat_dim == self.ndimension() - 1:
+            curr_idx = 0
+            res_list = []
+            index = [slice(None, None, None) for _ in range(rhs.ndimension())]
+            for t, size, rhs in zip(self.lazy_tensors, self.cat_dim_sizes, rhs_):
+                index[-2] = slice(curr_idx, curr_idx + size, None)
+                res_list.append(t._matmul(rhs[index]))
+                curr_idx += size
+            # copy result back to output device
+            res_list = [x.to(output_device) for x in res_list]
+            res = torch.sum(torch.stack(res_list), dim=0)
+        else:
+            while rhs.ndimension() < self.ndimension():
+                rhs = rhs.unsqueeze(0)
+            curr_idx = 0
+            res_list = []
+            index = [slice(None, None, None) for _ in range(self.ndimension())]
+            for t, size, rhs in zip(self.lazy_tensors, self.cat_dim_sizes, rhs_):
+                index[self.cat_dim] = slice(curr_idx, curr_idx + size, None)
+                res_list.append(t._matmul(rhs[index]))
+                curr_idx += size
+            # copy result back to output device
+            res_list = [x.to(output_device) for x in res_list]
+            res = torch.cat(res_list, dim=self.cat_dim)
+
+        if isvector:
+            res = res.squeeze(-1)
+        return res
+
+    def _size(self):
+        size = self.pre_cat_size + (self.cat_dim_len,) + self.post_cat_size
+        return torch.Size(size)
+
+    def _transpose_nonbatch(self):
+        if self.cat_dim == self.ndimension() - 2:
+            new_dim = self.cat_dim + 1
+        elif self.cat_dim == self.ndimension() - 1:
+            new_dim = self.cat_dim - 1
+        else:
+            new_dim = self.cat_dim
+        return self.__class__(*[t._transpose_nonbatch()
+                                for t in self.lazy_tensors], dim=new_dim, output_device=self.output_device)
+
+    def diag(self):
+        """
+        As :func:`torch.diag`, returns the diagonal of the matrix :math:`K` this LazyTensor represents as a vector.
+
+        Returns:
+            :obj:`torch.tensor`: The diagonal of :math:`K`. If :math:`K` is :math:`n \times n`, this will be a length
+            n vector. If this LazyTensor represents a batch (e.g., is :math:`b \times n \times n`), this will be a
+            :math:`b \times n` matrix of diagonals, one for each matrix in the batch.
+        """
+        size = self.size()
+        if size[-1] != size[-2]:
+            raise RuntimeError("Diag works on square matrices (or batches)")
+
+        row_col_iter = torch.arange(0, size[-1], dtype=torch.long)
+        if self.ndimension() == 3:
+            batch_iter = torch.arange(0, size[0], dtype=torch.long)
+            batch_iter = batch_iter.unsqueeze(1).repeat(1, size[1]).view(-1)
+            row_col_iter = row_col_iter.unsqueeze(1).repeat(size[0], 1).view(-1)
+            res = self._get_indices(row_col_iter, row_col_iter, batch_iter).view(size[0], size[1])
+        else:
+            res = self._get_indices(row_col_iter, row_col_iter)
+
+        return res.to(self.device)
+
+    def __getitem__(self, *indices):
+        res = super().__getitem__(*indices)
+        if not isinstance(res, CatLazyTensor):
+            res = res.to(self.device)
+
+        return res
+
+    def inv_matmul(self, right_tensor, left_tensor=None):
+        return super().inv_matmul(right_tensor, left_tensor).to(self.device)
+
+    def inv_quad(self, tensor):
+        return super().inv_quad(tensor).to(self.device)
+
+    def inv_quad_log_det(self, inv_quad_rhs=None, log_det=False, reduce_inv_quad=True):
+        res = super().inv_quad_log_det(inv_quad_rhs, log_det, reduce_inv_quad)
+        return tuple(r.to(self.device) for r in res)
+
+    def matmul(self, other):
+        return super().matmul(other).to(self.device)
+
+    @property
+    def device(self):
+        return self.output_device
+
+    @property
+    def devices(self):
+        return [t.device for t in self.lazy_tensors]
+
+    @property
+    def device_count(self):
+        return len(set(self.devices))
diff --git a/gpytorch/lazy/lazy_tensor.py b/gpytorch/lazy/lazy_tensor.py
index 1ec97796..af5dc066 100644
--- a/gpytorch/lazy/lazy_tensor.py
+++ b/gpytorch/lazy/lazy_tensor.py
@@ -67,7 +67,7 @@ class LazyTensor(object):
         operate on these batches if present.
     """
 
-    def _get_indices(self, *batch_indices, left_indices, right_indices):
+    def _get_indices(self, left_indices, right_indices, *batch_indices):
         """
         Returns entries of the matrix, indexed by batch, row, and column indices
         """
@@ -172,7 +172,7 @@ class LazyTensor(object):
             which does some additional work. Calling this method directly is discouraged.
 
         Args:
-            :attr:`indices` (tuple of int, slice, or LongTensor):
+            :attr:`indices` (tuple of `int`s, `slice`s, or `LongTensor`s):
                 A collection of indices for each of the dimensions. There will be exactly one index per dimension.
         """
         if settings.debug.on():


Test patch:
diff --git a/test/lazy/_lazy_tensor_test_case.py b/test/lazy/_lazy_tensor_test_case.py
index 5afbe145..a72c93a0 100644
--- a/test/lazy/_lazy_tensor_test_case.py
+++ b/test/lazy/_lazy_tensor_test_case.py
@@ -90,58 +90,74 @@ class RectangularLazyTensorTestCase(object):
         if lazy_tensor.ndimension() == 2:
             res = lazy_tensor[1]
             actual = evaluated[1]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             res = lazy_tensor[0:2].evaluate()
             actual = evaluated[0:2]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             res = lazy_tensor[:, 0:2].evaluate()
             actual = evaluated[:, 0:2]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             res = lazy_tensor[0:2, :].evaluate()
             actual = evaluated[0:2, :]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             res = lazy_tensor[..., 0:2].evaluate()
             actual = evaluated[..., 0:2]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             res = lazy_tensor[0:2, ...].evaluate()
             actual = evaluated[0:2, ...]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             res = lazy_tensor[..., 0:2, 2]
             actual = evaluated[..., 0:2, 2]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             res = lazy_tensor[0:2, ..., 2]
             actual = evaluated[0:2, ..., 2]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
 
         # Batch case
         else:
             res = lazy_tensor[1].evaluate()
             actual = evaluated[1]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             res = lazy_tensor[0:2].evaluate()
             actual = evaluated[0:2]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             res = lazy_tensor[:, 0:2].evaluate()
             actual = evaluated[:, 0:2]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
 
             for batch_index in product([1, slice(0, 2, None)], repeat=(lazy_tensor.dim() - 2)):
                 res = lazy_tensor.__getitem__((*batch_index, slice(0, 1, None), slice(0, 2, None))).evaluate()
                 actual = evaluated.__getitem__((*batch_index, slice(0, 1, None), slice(0, 2, None)))
+                self.assertEqual(res.shape, actual.shape)
                 self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
                 res = lazy_tensor.__getitem__((*batch_index, 1, slice(0, 2, None)))
                 actual = evaluated.__getitem__((*batch_index, 1, slice(0, 2, None)))
+                self.assertEqual(res.shape, actual.shape)
                 self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
                 res = lazy_tensor.__getitem__((*batch_index, slice(1, None, None), 2))
                 actual = evaluated.__getitem__((*batch_index, slice(1, None, None), 2))
+                self.assertEqual(res.shape, actual.shape)
                 self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
 
             # Ellipsis
             res = lazy_tensor.__getitem__((Ellipsis, slice(1, None, None), 2))
             actual = evaluated.__getitem__((Ellipsis, slice(1, None, None), 2))
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             res = lazy_tensor.__getitem__((slice(1, None, None), Ellipsis, 2))
             actual = evaluated.__getitem__((slice(1, None, None), Ellipsis, 2))
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
 
     def test_getitem_tensor_index(self):
@@ -152,21 +168,27 @@ class RectangularLazyTensorTestCase(object):
         if lazy_tensor.ndimension() == 2:
             index = (torch.tensor([0, 0, 1, 2]), torch.tensor([0, 1, 0, 2]))
             res, actual = lazy_tensor[index], evaluated[index]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             index = (torch.tensor([0, 0, 1, 2]), slice(None, None, None))
-            res, actual = lazy_tensor[index], evaluated[index]
+            res, actual = gpytorch.delazify(lazy_tensor[index]), evaluated[index]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             index = (slice(None, None, None), torch.tensor([0, 0, 1, 2]))
-            res, actual = lazy_tensor[index], evaluated[index]
+            res, actual = gpytorch.delazify(lazy_tensor[index]), evaluated[index]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             index = (torch.tensor([0, 0, 1, 2]), Ellipsis)
-            res, actual = lazy_tensor[index], evaluated[index]
+            res, actual = gpytorch.delazify(lazy_tensor[index]), evaluated[index]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             index = (Ellipsis, torch.tensor([0, 0, 1, 2]))
-            res, actual = lazy_tensor[index], evaluated[index]
+            res, actual = gpytorch.delazify(lazy_tensor[index]), evaluated[index]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
             index = (Ellipsis, torch.tensor([0, 0, 1, 2]), torch.tensor([0, 1, 0, 2]))
             res, actual = lazy_tensor[index], evaluated[index]
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
 
         # Batch case
@@ -176,23 +198,31 @@ class RectangularLazyTensorTestCase(object):
             ):
                 index = (*batch_index, torch.tensor([0, 1, 0, 2]), torch.tensor([1, 2, 0, 1]))
                 res, actual = lazy_tensor[index], evaluated[index]
+                self.assertEqual(res.shape, actual.shape)
                 self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
                 index = (*batch_index, torch.tensor([0, 1, 0, 2]), slice(None, None, None))
-                res, actual = lazy_tensor[index], evaluated[index]
+                res, actual = gpytorch.delazify(lazy_tensor[index]), evaluated[index]
+                self.assertEqual(res.shape, actual.shape)
                 self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
                 index = (*batch_index, slice(None, None, None), torch.tensor([0, 1, 2, 1]))
-                res, actual = lazy_tensor[index], evaluated[index]
+                res, actual = gpytorch.delazify(lazy_tensor[index]), evaluated[index]
+                self.assertEqual(res.shape, actual.shape)
                 self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
                 index = (*batch_index, slice(None, None, None), slice(None, None, None))
                 res, actual = lazy_tensor[index].evaluate(), evaluated[index]
+                self.assertEqual(res.shape, actual.shape)
                 self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
 
             # Ellipsis
             res = lazy_tensor.__getitem__((Ellipsis, torch.tensor([0, 1, 0, 2]), torch.tensor([1, 2, 0, 1])))
             actual = evaluated.__getitem__((Ellipsis, torch.tensor([0, 1, 0, 2]), torch.tensor([1, 2, 0, 1])))
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
-            res = lazy_tensor.__getitem__((torch.tensor([0, 1, 0, 1]), Ellipsis, torch.tensor([1, 2, 0, 1])))
+            res = gpytorch.delazify(
+                lazy_tensor.__getitem__((torch.tensor([0, 1, 0, 1]), Ellipsis, torch.tensor([1, 2, 0, 1])))
+            )
             actual = evaluated.__getitem__((torch.tensor([0, 1, 0, 1]), Ellipsis, torch.tensor([1, 2, 0, 1])))
+            self.assertEqual(res.shape, actual.shape)
             self.assertLess(((res - actual).abs() / actual.abs().clamp(1, 1e5)).max().item(), 1e-1)
 
 
diff --git a/test/lazy/test_cat_lazy_tensor.py b/test/lazy/test_cat_lazy_tensor.py
new file mode 100644
index 00000000..e079a5a6
--- /dev/null
+++ b/test/lazy/test_cat_lazy_tensor.py
@@ -0,0 +1,54 @@
+#!/usr/bin/env python3
+
+import torch
+import unittest
+from gpytorch.lazy import CatLazyTensor, NonLazyTensor
+from test.lazy._lazy_tensor_test_case import LazyTensorTestCase
+
+
+class TestCatLazyTensor(LazyTensorTestCase, unittest.TestCase):
+    seed = 0
+    should_test_sample = True
+
+    def create_lazy_tensor(self):
+        root = torch.randn(6, 7)
+        self.psd_mat = root.matmul(root.t())
+
+        slice1_mat = self.psd_mat[:2, :].requires_grad_()
+        slice2_mat = self.psd_mat[2:4, :].requires_grad_()
+        slice3_mat = self.psd_mat[4:6, :].requires_grad_()
+
+        slice1 = NonLazyTensor(slice1_mat)
+        slice2 = NonLazyTensor(slice2_mat)
+        slice3 = NonLazyTensor(slice3_mat)
+
+        return CatLazyTensor(slice1, slice2, slice3, dim=-2)
+
+    def evaluate_lazy_tensor(self, lazy_tensor):
+        return self.psd_mat.detach().clone().requires_grad_()
+
+
+class TestCatLazyTensorBatch(LazyTensorTestCase, unittest.TestCase):
+    seed = 0
+    should_test_sample = True
+
+    def create_lazy_tensor(self):
+        root = torch.randn(3, 6, 7)
+        self.psd_mat = root.matmul(root.transpose(-2, -1))
+
+        slice1_mat = self.psd_mat[..., :2, :].requires_grad_()
+        slice2_mat = self.psd_mat[..., 2:4, :].requires_grad_()
+        slice3_mat = self.psd_mat[..., 4:6, :].requires_grad_()
+
+        slice1 = NonLazyTensor(slice1_mat)
+        slice2 = NonLazyTensor(slice2_mat)
+        slice3 = NonLazyTensor(slice3_mat)
+
+        return CatLazyTensor(slice1, slice2, slice3, dim=-2)
+
+    def evaluate_lazy_tensor(self, lazy_tensor):
+        return self.psd_mat.detach().clone().requires_grad_()
+
+
+if __name__ == "__main__":
+    unittest.main()