Repository: cornellius-gp/gpytorch
Pull request #1057 (2020-02-24 22:06:24)

Orthogonally decoupled variational strategy
Implements [Orthogonally Decoupled Variational Gaussian Processes](https://arxiv.org/abs/1809.08820) from Salimbeni et al. (2018).

jacobrgardner: Actual test failures are obviously unrelated to this PR -- probably something changed on pytorch master?

We should also fix our `add_` deprecation warnings before cutting our next minor release.
gpleiss: @jacobrgardner - yeah I'm fixing the test failures right now.
gpleiss: @jacobrgardner #1058 fixes the error from pytorch#master

Patch:
diff --git a/docs/source/variational.rst b/docs/source/variational.rst
index ec961876..78e7227a 100644
--- a/docs/source/variational.rst
+++ b/docs/source/variational.rst
@@ -63,6 +63,13 @@ deterministic relationship between :math:`\mathbf f` and :math:`\mathbf u`.
    :members:
 
 
+:hidden:`OrthogonallyDecoupledVariationalStrategy`
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. autoclass:: OrthogonallyDecoupledVariationalStrategy
+   :members:
+
+
 :hidden:`UnwhitenedVariationalStrategy`
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/examples/04_Variational_and_Approximate_GPs/Modifying_the_variational_strategy_and_distribution.ipynb b/examples/04_Variational_and_Approximate_GPs/Modifying_the_variational_strategy_and_distribution.ipynb
new file mode 100644
index 00000000..b6016fa7
--- /dev/null
+++ b/examples/04_Variational_and_Approximate_GPs/Modifying_the_variational_strategy_and_distribution.ipynb
@@ -0,0 +1,476 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import math\n",
+    "import torch\n",
+    "import gpytorch\n",
+    "import tqdm\n",
+    "import matplotlib.pyplot as plt\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Modifying the Variational Strategy/Variational Distribution\n",
+    "\n",
+    "The predictive distribution for approximate GPs is given by\n",
+    "\n",
+    "$$\n",
+    "p( \\mathbf f(\\mathbf x^*) ) = \\int_{\\mathbf u} p( f(\\mathbf x^*) \\mid \\mathbf u) \\: q(\\mathbf u) \\: d\\mathbf u,\n",
+    "\\quad\n",
+    "q(\\mathbf u) = \\mathcal N( \\mathbf m, \\mathbf S).\n",
+    "$$\n",
+    "\n",
+    "$\\mathbf u$ represents the function values at the $m$ inducing points.\n",
+    "Here, $\\mathbf m \\in \\mathbb R^m$ and $\\mathbf S \\in \\mathbb R^{m \\times m}$ are learnable parameters.\n",
+    "\n",
+    "If $m$ (the number of inducing points) is quite large, the number of learnable parameters in $\\mathbf S$ can be quite unwieldy.\n",
+    "Furthermore, a large $m$ might make some of the computations rather slow.\n",
+    "Here we show a few ways to use different [variational distributions](https://gpytorch.readthedocs.io/en/latest/variational.html#variational-distributions) and\n",
+    "[variational strategies](https://gpytorch.readthedocs.io/en/latest/variational.html#variational-strategies) to accomplish this."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Experimental setup\n",
+    "\n",
+    "We're going to train an approximate GP on a medium-sized regression dataset, taken from the UCI repository."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import urllib.request\n",
+    "import os\n",
+    "from scipy.io import loadmat\n",
+    "from math import floor\n",
+    "\n",
+    "\n",
+    "# this is for running the notebook in our testing framework\n",
+    "smoke_test = ('CI' in os.environ)\n",
+    "\n",
+    "\n",
+    "if not smoke_test and not os.path.isfile('../elevators.mat'):\n",
+    "    print('Downloading \\'elevators\\' UCI dataset...')\n",
+    "    urllib.request.urlretrieve('https://drive.google.com/uc?export=download&id=1jhWL3YUHvXIaftia4qeAyDwVxo6j1alk', '../elevators.mat')\n",
+    "\n",
+    "\n",
+    "if smoke_test:  # this is for running the notebook in our testing framework\n",
+    "    X, y = torch.randn(1000, 3), torch.randn(1000)\n",
+    "else:\n",
+    "    data = torch.Tensor(loadmat('../elevators.mat')['data'])\n",
+    "    X = data[:, :-1]\n",
+    "    X = X - X.min(0)[0]\n",
+    "    X = 2 * (X / X.max(0)[0]) - 1\n",
+    "    y = data[:, -1]\n",
+    "\n",
+    "\n",
+    "train_n = int(floor(0.8 * len(X)))\n",
+    "train_x = X[:train_n, :].contiguous()\n",
+    "train_y = y[:train_n].contiguous()\n",
+    "\n",
+    "test_x = X[train_n:, :].contiguous()\n",
+    "test_y = y[train_n:].contiguous()\n",
+    "\n",
+    "if torch.cuda.is_available():\n",
+    "    train_x, train_y, test_x, test_y = train_x.cuda(), train_y.cuda(), test_x.cuda(), test_y.cuda()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from torch.utils.data import TensorDataset, DataLoader\n",
+    "train_dataset = TensorDataset(train_x, train_y)\n",
+    "train_loader = DataLoader(train_dataset, batch_size=500, shuffle=True)\n",
+    "\n",
+    "test_dataset = TensorDataset(test_x, test_y)\n",
+    "test_loader = DataLoader(test_dataset, batch_size=500, shuffle=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Some quick training/testing code\n",
+    "\n",
+    "This will allow us to train/test different model classes."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# this is for running the notebook in our testing framework\n",
+    "num_epochs = 1 if smoke_test else 10\n",
+    "\n",
+    "\n",
+    "# Our testing script takes in a GPyTorch MLL (objective function) class\n",
+    "# and then trains/tests an approximate GP with it on the supplied dataset\n",
+    "\n",
+    "def train_and_test_approximate_gp(model_cls):\n",
+    "    inducing_points = torch.randn(128, train_x.size(-1), dtype=train_x.dtype, device=train_x.device)\n",
+    "    model = model_cls(inducing_points)\n",
+    "    likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
+    "    mll = gpytorch.mlls.VariationalELBO(likelihood, model, num_data=train_y.numel())\n",
+    "    optimizer = torch.optim.Adam(list(model.parameters()) + list(likelihood.parameters()), lr=0.1)\n",
+    "    \n",
+    "    if torch.cuda.is_available():\n",
+    "        model = model.cuda()\n",
+    "        likelihood = likelihood.cuda()\n",
+    "\n",
+    "    # Training\n",
+    "    model.train()\n",
+    "    likelihood.train()\n",
+    "    epochs_iter = tqdm.notebook.tqdm(range(num_epochs), desc=f\"Training {model_cls.__name__}\")\n",
+    "    for i in epochs_iter:\n",
+    "        # Within each iteration, we will go over each minibatch of data\n",
+    "        for x_batch, y_batch in train_loader:\n",
+    "            optimizer.zero_grad()\n",
+    "            output = model(x_batch)\n",
+    "            loss = -mll(output, y_batch)\n",
+    "            epochs_iter.set_postfix(loss=loss.item())\n",
+    "            loss.backward()\n",
+    "            optimizer.step()\n",
+    "            \n",
+    "    # Testing\n",
+    "    model.eval()\n",
+    "    likelihood.eval()\n",
+    "    means = torch.tensor([0.])\n",
+    "    with torch.no_grad():\n",
+    "        for x_batch, y_batch in test_loader:\n",
+    "            preds = model(x_batch)\n",
+    "            means = torch.cat([means, preds.mean.cpu()])\n",
+    "    means = means[1:]\n",
+    "    error = torch.mean(torch.abs(means - test_y.cpu()))\n",
+    "    print(f\"Test {model_cls.__name__} MAE: {error.item()}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The Standard Approach\n",
+    "\n",
+    "As a default, we'll use the default [VariationalStrategy](https://gpytorch.readthedocs.io/en/latest/variational.html#id1) class with a [CholeskyVariationalDistribution](https://gpytorch.readthedocs.io/en/latest/variational.html#choleskyvariationaldistribution).\n",
+    "The `CholeskyVariationalDistribution` class allows $\\mathbf S$ to be on any positive semidefinite matrix. This is the most general/expressive option for approximate GPs."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class StandardApproximateGP(gpytorch.models.ApproximateGP):\n",
+    "    def __init__(self, inducing_points):\n",
+    "        variational_distribution = gpytorch.variational.CholeskyVariationalDistribution(inducing_points.size(-2))\n",
+    "        variational_strategy = gpytorch.variational.VariationalStrategy(\n",
+    "            self, inducing_points, variational_distribution, learn_inducing_locations=True\n",
+    "        )\n",
+    "        super().__init__(variational_strategy)\n",
+    "        self.mean_module = gpytorch.means.ConstantMean()\n",
+    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\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)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "7ecb0000e03d4b95b5e9dd762da0c1bb",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "HBox(children=(IntProgress(value=0, description='Training StandardApproximateGP', max=10, style=ProgressStyle(…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "Test StandardApproximateGP MAE: 0.10098349303007126\n"
+     ]
+    }
+   ],
+   "source": [
+    "train_and_test_approximate_gp(StandardApproximateGP)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Reducing parameters\n",
+    "\n",
+    "### MeanFieldVariationalDistribution: a diagonal $\\mathbf S$ matrix \n",
+    "\n",
+    "One way to reduce the number of parameters is to restrict that $\\mathbf S$ is only diagonal. This is less expressive, but the number of parameters is now linear in $m$ instead of quadratic.\n",
+    "\n",
+    "All we have to do is take the previous example, and change `CholeskyVariationalDistribution` (full $\\mathbf S$ matrix) to [MeanFieldVariationalDistribution](https://gpytorch.readthedocs.io/en/latest/variational.html#meanfieldvariationaldistribution) (diagonal $\\mathbf S$ matrix)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class MeanFieldApproximateGP(gpytorch.models.ApproximateGP):\n",
+    "    def __init__(self, inducing_points):\n",
+    "        variational_distribution = gpytorch.variational.MeanFieldVariationalDistribution(inducing_points.size(-2))\n",
+    "        variational_strategy = gpytorch.variational.VariationalStrategy(\n",
+    "            self, inducing_points, variational_distribution, learn_inducing_locations=True\n",
+    "        )\n",
+    "        super().__init__(variational_strategy)\n",
+    "        self.mean_module = gpytorch.means.ConstantMean()\n",
+    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\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)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 21,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "2e95f8db268647f4b069758d9ed9dffd",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "HBox(children=(IntProgress(value=0, description='Training MeanFieldApproximateGP', max=10, style=ProgressStyle…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "Test MeanFieldApproximateGP MAE: 0.07848489284515381\n"
+     ]
+    }
+   ],
+   "source": [
+    "train_and_test_approximate_gp(MeanFieldApproximateGP)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### DeltaVariationalDistribution: no $\\mathbf S$ matrix \n",
+    "\n",
+    "A more extreme method of reducing parameters is to get rid of $\\mathbf S$ entirely. This corresponds to learning a delta distribution ($\\mathbf u = \\mathbf m$) rather than a multivariate Normal distribution for $\\mathbf u$. In other words, this corresponds to performing MAP estimation rather than variational inference.\n",
+    "\n",
+    "In GPyTorch, getting rid of $\\mathbf S$ can be accomplished by using a [DeltaVariationalDistribution](https://gpytorch.readthedocs.io/en/latest/variational.html#deltavariationaldistribution)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class MAPApproximateGP(gpytorch.models.ApproximateGP):\n",
+    "    def __init__(self, inducing_points):\n",
+    "        variational_distribution = gpytorch.variational.DeltaVariationalDistribution(inducing_points.size(-2))\n",
+    "        variational_strategy = gpytorch.variational.VariationalStrategy(\n",
+    "            self, inducing_points, variational_distribution, learn_inducing_locations=True\n",
+    "        )\n",
+    "        super().__init__(variational_strategy)\n",
+    "        self.mean_module = gpytorch.means.ConstantMean()\n",
+    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\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)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "b6cb9a96daea4abc8d092f82d38c0c40",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "HBox(children=(IntProgress(value=0, description='Training MAPApproximateGP', max=10, style=ProgressStyle(descr…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "Test MAPApproximateGP MAE: 0.08846496045589447\n"
+     ]
+    }
+   ],
+   "source": [
+    "train_and_test_approximate_gp(MAPApproximateGP)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Reducing computation (through decoupled inducing points)\n",
+    "\n",
+    "One way to reduce the computational complexity is to use separate inducing points for the mean and covariance computations. The [Orthogonally Decoupled Variational Gaussian Processes](https://arxiv.org/abs/1809.08820) method of Salimbeni et al. (2018) uses more inducing points for the (computationally easy) mean computations and fewer inducing points for the (computationally intensive) covariance computations.\n",
+    "\n",
+    "In GPyTorch we implement this method in a modular way. The [OrthogonallyDecoupledVariationalStrategy](http://gpytorch.ai/variational.html#gpytorch.variational.OrthogonallyDecoupledVariationalStrategy) defines the variational strategy for the mean inducing points. It wraps an existing variational strategy/distribution that defines the covariance inducing points:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 28,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def make_orthogonal_vs(model, train_x):\n",
+    "    mean_inducing_points = torch.randn(1000, train_x.size(-1), dtype=train_x.dtype, device=train_x.device)\n",
+    "    covar_inducing_points = torch.randn(100, train_x.size(-1), dtype=train_x.dtype, device=train_x.device)\n",
+    "\n",
+    "    covar_variational_strategy = gpytorch.variational.VariationalStrategy(\n",
+    "        model, covar_inducing_points,\n",
+    "        gpytorch.variational.CholeskyVariationalDistribution(covar_inducing_points.size(-2)),\n",
+    "        learn_inducing_locations=True\n",
+    "    )\n",
+    "\n",
+    "    variational_strategy = gpytorch.variational.OrthogonallyDecoupledVariationalStrategy(\n",
+    "        covar_variational_strategy, mean_inducing_points,\n",
+    "        gpytorch.variational.DeltaVariationalDistribution(mean_inducing_points.size(-2)),\n",
+    "    )\n",
+    "    return variational_strategy"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Putting it all together we have:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 29,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class OrthDecoupledApproximateGP(gpytorch.models.ApproximateGP):\n",
+    "    def __init__(self, inducing_points):\n",
+    "        variational_distribution = gpytorch.variational.DeltaVariationalDistribution(inducing_points.size(-2))\n",
+    "        variational_strategy = make_orthogonal_vs(self, train_x)\n",
+    "        super().__init__(variational_strategy)\n",
+    "        self.mean_module = gpytorch.means.ConstantMean()\n",
+    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\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)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 30,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "ef883369bfc84b919c2999823052f04b",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "HBox(children=(IntProgress(value=0, description='Training OrthDecoupledApproximateGP', max=10, style=ProgressS…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "Test OrthDecoupledApproximateGP MAE: 0.08162340521812439\n"
+     ]
+    }
+   ],
+   "source": [
+    "train_and_test_approximate_gp(OrthDecoupledApproximateGP)"
+   ]
+  }
+ ],
+ "metadata": {
+  "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.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/examples/04_Variational_and_Approximate_GPs/index.rst b/examples/04_Variational_and_Approximate_GPs/index.rst
index b6410c8b..702f488a 100644
--- a/examples/04_Variational_and_Approximate_GPs/index.rst
+++ b/examples/04_Variational_and_Approximate_GPs/index.rst
@@ -18,10 +18,12 @@ All approximate models consist of the following 3 composible objects:
 - :obj:`~gpytorch.mlls._ApproximateMarginalLogLikelihood`, which defines the objective function
   to learn the approximate posterior (e.g. variational ELBO).
 
+(See the `strategy/distribution comparison`_ for examples of the different classes.)
 The variational documentation has more information on how to use these objects.
 Here we provide some examples which highlight some of the common use cases:
 
 - **Large-scale regression** (when exact methods are too memory intensive): see the `stochastic variational regression example`_.
+- **Variational distribution options** for different scalability/expressiveness: see the `strategy/distribution comparison`_.
 - **Alternative optimization objectives** for the GP's predictive distribution: see the `approximate GP objective functions notebook`_.
   This example compares and contrasts the variational ELBO with the predictive log likelihood of Jankowiak et al., 2019.
 - **Classification**: see the `non-Gaussian likelihood notebook`_.
@@ -33,11 +35,15 @@ Here we provide some examples which highlight some of the common use cases:
    :hidden:
 
    SVGP_Regression_CUDA.ipynb
+   Modifying_the_variational_strategy_and_distribution.ipynb
    Approximate_GP_Objective_Functions.ipynb
    Non_Gaussian_Likelihoods.ipynb
    SVGP_Multitask_GP_Regression.ipynb
    GP_Regression_with_Uncertain_Inputs.ipynb
 
+.. _strategy/distribution comparison:
+  ./Modifying_the_variational_strategy_and_distribution.ipynb
+
 .. _stochastic variational regression example:
   ./SVGP_Regression_CUDA.ipynb
 
diff --git a/gpytorch/variational/__init__.py b/gpytorch/variational/__init__.py
index 01c27ecf..4177a348 100644
--- a/gpytorch/variational/__init__.py
+++ b/gpytorch/variational/__init__.py
@@ -8,6 +8,7 @@ from .delta_variational_distribution import DeltaVariationalDistribution
 from .grid_interpolation_variational_strategy import GridInterpolationVariationalStrategy
 from .mean_field_variational_distribution import MeanFieldVariationalDistribution
 from .multitask_variational_strategy import MultitaskVariationalStrategy
+from .orthogonally_decoupled_variational_strategy import OrthogonallyDecoupledVariationalStrategy
 from .unwhitened_variational_strategy import UnwhitenedVariationalStrategy
 from .variational_strategy import VariationalStrategy
 from .whitened_variational_strategy import WhitenedVariationalStrategy
@@ -17,6 +18,7 @@ __all__ = [
     "AdditiveGridInterpolationVariationalStrategy",
     "GridInterpolationVariationalStrategy",
     "MultitaskVariationalStrategy",
+    "OrthogonallyDecoupledVariationalStrategy",
     "VariationalStrategy",
     "UnwhitenedVariationalStrategy",
     "WhitenedVariationalStrategy",
diff --git a/gpytorch/variational/orthogonally_decoupled_variational_strategy.py b/gpytorch/variational/orthogonally_decoupled_variational_strategy.py
new file mode 100644
index 00000000..95457055
--- /dev/null
+++ b/gpytorch/variational/orthogonally_decoupled_variational_strategy.py
@@ -0,0 +1,86 @@
+#!/usr/bin/env python3
+
+import torch
+
+from ..distributions import MultivariateNormal
+from ..utils.memoize import cached
+from ._variational_strategy import _VariationalStrategy
+from .delta_variational_distribution import DeltaVariationalDistribution
+
+
+class OrthogonallyDecoupledVariationalStrategy(_VariationalStrategy):
+    r"""
+    Implements orthogonally decoupled VGPs as defined in `Salimbeni et al. (2018)`_.
+    This variational strategy uses a different set of inducing points for the mean and covariance functions.
+    The idea is to use more inducing points for the (computationally efficient) mean and fewer inducing points for the
+    (computationally expensive) covaraince.
+
+    This variational strategy defines the inducing points/:obj:`~gpytorch.variational._VariationalDistribution`
+    for the mean function.
+    It then wraps a different :obj:`~gpytorch.variational._VariationalStrategy` which
+    defines the covariance inducing points.
+
+    Example:
+        >>> mean_inducing_points = torch.randn(1000, train_x.size(-1), dtype=train_x.dtype, device=train_x.device)
+        >>> covar_inducing_points = torch.randn(100, train_x.size(-1), dtype=train_x.dtype, device=train_x.device)
+        >>>
+        >>> covar_variational_strategy = gpytorch.variational.VariationalStrategy(
+        >>>     model, covar_inducing_points,
+        >>>     gpytorch.variational.CholeskyVariationalDistribution(covar_inducing_points.size(-2)),
+        >>>     learn_inducing_locations=True
+        >>> )
+        >>>
+        >>> variational_strategy = gpytorch.variational.OrthogonallyDecoupledVariationalStrategy(
+        >>>     covar_variational_strategy, mean_inducing_points,
+        >>>     gpytorch.variational.DeltaVariationalDistribution(mean_inducing_points.size(-2)),
+        >>> )
+
+    .. _Salimbeni et al. (2018):
+        https://arxiv.org/abs/1809.08820
+    """
+
+    def __init__(self, model, inducing_points, variational_distribution):
+        if not isinstance(variational_distribution, DeltaVariationalDistribution):
+            raise NotImplementedError(
+                "OrthogonallyDecoupledVariationalStrategy currently works with DeltaVariationalDistribution"
+            )
+
+        super().__init__(model, inducing_points, variational_distribution, learn_inducing_locations=True)
+        self.base_variational_strategy = model
+
+    @property
+    @cached(name="prior_distribution_memo")
+    def prior_distribution(self):
+        out = self.model(self.inducing_points)
+        res = MultivariateNormal(out.mean, out.lazy_covariance_matrix.add_jitter())
+        return res
+
+    def forward(self, x, inducing_points, inducing_values, variational_inducing_covar=None):
+        if variational_inducing_covar is not None:
+            raise NotImplementedError(
+                "OrthogonallyDecoupledVariationalStrategy currently works with DeltaVariationalDistribution"
+            )
+
+        num_data = x.size(-2)
+        full_output = self.model(torch.cat([x, inducing_points], dim=-2))
+        full_mean = full_output.mean
+        full_covar = full_output.lazy_covariance_matrix
+
+        if self.training:
+            induc_mean = full_mean[..., num_data:]
+            induc_induc_covar = full_covar[..., num_data:, num_data:]
+            self._memoize_cache["prior_distribution_memo"] = MultivariateNormal(induc_mean, induc_induc_covar)
+
+        test_mean = full_mean[..., :num_data]
+        data_induc_covar = full_covar[..., :num_data, num_data:]
+        predictive_mean = (data_induc_covar @ inducing_values.unsqueeze(-1)).squeeze(-1).add(test_mean)
+        predictive_covar = full_covar[..., :num_data, :num_data]
+
+        # Return the distribution
+        return MultivariateNormal(predictive_mean, predictive_covar)
+
+    def kl_divergence(self):
+        mean = self.variational_distribution.mean
+        induc_induc_covar = self.prior_distribution.lazy_covariance_matrix
+        kl = self.model.kl_divergence() + ((induc_induc_covar @ mean.unsqueeze(-1)).squeeze(-1) * mean).sum(-1).mul(0.5)
+        return kl


Test patch:
diff --git a/test/examples/test_decoupled_svgp_regression.py b/test/examples/test_decoupled_svgp_regression.py
new file mode 100644
index 00000000..0bfd2f83
--- /dev/null
+++ b/test/examples/test_decoupled_svgp_regression.py
@@ -0,0 +1,99 @@
+#!/usr/bin/env python3
+
+import math
+import warnings
+import unittest
+from unittest.mock import MagicMock, patch
+
+import gpytorch
+import torch
+from gpytorch.likelihoods import GaussianLikelihood
+from gpytorch.models import ApproximateGP
+from gpytorch.test.base_test_case import BaseTestCase
+from gpytorch.lazy import ExtraComputationWarning
+from torch import optim
+
+
+def train_data():
+    train_x = torch.linspace(0, 1, 260)
+    train_y = torch.cos(train_x * (2 * math.pi))
+    return train_x, train_y
+
+
+class SVGPRegressionModel(ApproximateGP):
+    def __init__(self, inducing_points, base_inducing_points):
+        base_variational_distribution = gpytorch.variational.CholeskyVariationalDistribution(
+            base_inducing_points.size(-1)
+        )
+        variational_distribution = gpytorch.variational.DeltaVariationalDistribution(inducing_points.size(-1))
+        variational_strategy = gpytorch.variational.OrthogonallyDecoupledVariationalStrategy(
+            gpytorch.variational.VariationalStrategy(
+                self, base_inducing_points, base_variational_distribution, learn_inducing_locations=True,
+            ), inducing_points, variational_distribution
+        )
+        super(SVGPRegressionModel, self).__init__(variational_strategy)
+        self.mean_module = gpytorch.means.ConstantMean()
+        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())
+
+    def forward(self, x):
+        mean_x = self.mean_module(x)
+        covar_x = self.covar_module(x)
+        latent_pred = gpytorch.distributions.MultivariateNormal(mean_x, covar_x)
+        return latent_pred
+
+
+class TestSVGPRegression(BaseTestCase, unittest.TestCase):
+    seed = 0
+
+    def test_regression_error(
+        self,
+        mll_cls=gpytorch.mlls.VariationalELBO,
+        distribution_cls=gpytorch.variational.CholeskyVariationalDistribution,
+    ):
+        train_x, train_y = train_data()
+        likelihood = GaussianLikelihood()
+        model = SVGPRegressionModel(torch.linspace(0, 1, 128), torch.linspace(0, 1, 16))
+        mll = mll_cls(likelihood, model, num_data=len(train_y))
+
+        # Find optimal model hyperparameters
+        model.train()
+        likelihood.train()
+        optimizer = optim.Adam([{"params": model.parameters()}, {"params": likelihood.parameters()}], lr=0.01)
+
+        _wrapped_cg = MagicMock(wraps=gpytorch.utils.linear_cg)
+        _cg_mock = patch("gpytorch.utils.linear_cg", new=_wrapped_cg)
+        with warnings.catch_warnings(record=True) as ws, _cg_mock as cg_mock:
+            for _ in range(75):
+                optimizer.zero_grad()
+                output = model(train_x)
+                loss = -mll(output, train_y)
+                loss.backward()
+                optimizer.step()
+
+            for param in model.parameters():
+                self.assertTrue(param.grad is not None)
+                self.assertGreater(param.grad.norm().item(), 0)
+            for param in likelihood.parameters():
+                self.assertTrue(param.grad is not None)
+                self.assertGreater(param.grad.norm().item(), 0)
+
+            # Set back to eval mode
+            model.eval()
+            likelihood.eval()
+            test_preds = likelihood(model(train_x)).mean.squeeze()
+            mean_abs_error = torch.mean(torch.abs(train_y - test_preds) / 2)
+            self.assertLess(mean_abs_error.item(), 1e-1)
+
+            # Make sure CG was called (or not), and no warnings were thrown
+            self.assertFalse(cg_mock.called)
+            self.assertFalse(any(issubclass(w.category, ExtraComputationWarning) for w in ws))
+
+    def test_predictive_ll_regression_error(self):
+        return self.test_regression_error(
+            mll_cls=gpytorch.mlls.PredictiveLogLikelihood,
+            distribution_cls=gpytorch.variational.MeanFieldVariationalDistribution,
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/variational/test_orthogonally_decoupled_variational_strategy.py b/test/variational/test_orthogonally_decoupled_variational_strategy.py
new file mode 100644
index 00000000..0d84d436
--- /dev/null
+++ b/test/variational/test_orthogonally_decoupled_variational_strategy.py
@@ -0,0 +1,72 @@
+#!/usr/bin/env python3
+
+import unittest
+
+import torch
+
+import gpytorch
+from gpytorch.test.variational_test_case import VariationalTestCase
+
+
+def likelihood_cls():
+    return gpytorch.likelihoods.GaussianLikelihood()
+
+
+def strategy_cls(model, inducing_points, variational_distribution, learn_inducing_locations):
+    base_inducing_points = torch.randn(8, inducing_points.size(-1), device=inducing_points.device)
+    base_variational_distribution = gpytorch.variational.CholeskyVariationalDistribution(8)
+    return gpytorch.variational.OrthogonallyDecoupledVariationalStrategy(
+        gpytorch.variational.VariationalStrategy(
+            model, base_inducing_points, base_variational_distribution, learn_inducing_locations
+        ),
+        inducing_points,
+        variational_distribution,
+    )
+
+
+class TestOrthogonallyDecoupledVariationalGP(VariationalTestCase, unittest.TestCase):
+    @property
+    def batch_shape(self):
+        return torch.Size([])
+
+    @property
+    def distribution_cls(self):
+        return gpytorch.variational.DeltaVariationalDistribution
+
+    @property
+    def likelihood_cls(self):
+        return likelihood_cls
+
+    @property
+    def mll_cls(self):
+        return gpytorch.mlls.VariationalELBO
+
+    @property
+    def strategy_cls(self):
+        return strategy_cls
+
+    def test_training_iteration(self, *args, **kwargs):
+        cg_mock, cholesky_mock = super().test_training_iteration(*args, **kwargs)
+        self.assertFalse(cg_mock.called)
+        self.assertEqual(cholesky_mock.call_count, 3)  # One for each forward pass, and for computing prior dist
+
+    def test_eval_iteration(self, *args, **kwargs):
+        cg_mock, cholesky_mock = super().test_eval_iteration(*args, **kwargs)
+        self.assertFalse(cg_mock.called)
+        self.assertEqual(cholesky_mock.call_count, 1)  # One to compute cache, that's it!
+
+
+class TestOrthogonallyDecoupledPredictiveGP(TestOrthogonallyDecoupledVariationalGP):
+    @property
+    def mll_cls(self):
+        return gpytorch.mlls.PredictiveLogLikelihood
+
+
+class TestOrthogonallyDecoupledRobustVGP(TestOrthogonallyDecoupledVariationalGP):
+    @property
+    def mll_cls(self):
+        return gpytorch.mlls.GammaRobustVariationalELBO
+
+
+if __name__ == "__main__":
+    unittest.main()