We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
[Docs] Deep Sigma Point Processes
Hi,

The DSPPs documentaiton makes it seem like theyare coded up very similarly to the standard Deep GPs, and this does seem true for[ two layers](https://docs.gpytorch.ai/en/stable/examples/05_Deep_Gaussian_Processes/Deep_Sigma_Point_Processes.html). 

However, if I try to add a third layer (a second hidden layer) and update the input and output dims and make sure the forward call is correct, the loss function returns a vector rather than a single value (which the optimsier obviously doesn't like).

I'm guessing these losses are meant to be weighted according to the quadrature points. Is this true? Perhpas this could be more clear in the documentation (or alternatively, should this be done automatically in the backend?)

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `gpytorch/models/deep_gps/dspp.py`
Content:
```
1 import torch
2 
3 from gpytorch import settings
4 from gpytorch.distributions import MultitaskMultivariateNormal
5 from gpytorch.lazy import BlockDiagLazyTensor
6 
7 from ..approximate_gp import ApproximateGP
8 from .deep_gp import DeepGP, DeepGPLayer
9 
10 
11 class DSPPLayer(DeepGPLayer):
12     """
13     Represents a layer in a DSPP where inference is performed using the techniques defined in Jankowiak et al., 2020.
14 
15     From an end user perspective, the functionality and usage of a DSPPLayer is essentially identical to that of a
16     DeepGPLayer. It is therefore recommended that you review the documentation for DeepGPLayer.
17 
18     :param ~gpytorch.variational.VariationalStrategy variational_strategy: Strategy for
19         changing q(u) -> q(f) (see other VI docs)
20     :param int input_dims: Dimensionality of input data expected by each GP
21     :param int output_dims: (default None) Number of GPs in this layer, equivalent to
22         output dimensionality. If set to `None`, then the output dimension will be squashed.
23     :param int num_quad_sites: Number of quadrature sites to use. Also the number of Gaussians in the mixture output
24         by this layer.
25 
26     Again, refer to the documentation for DeepGPLayer or our example notebooks for full details on what calling a
27     DSPPLayer module does. The high level overview is that if a tensor `x` is `n x d` then
28 
29         >>> hidden_gp2(hidden_gp1(x))
30 
31     will return a `num_quad_sites` by `output_dims` set of Gaussians, where for each output dim the first batch dim
32     represents a weighted mixture of `num_quad_sites` Gaussians with weights given by DSPP.quad_weights (see DSPP below)
33     """
34 
35     def __init__(self, variational_strategy, input_dims, output_dims, num_quad_sites=3, quad_sites=None):
36         super().__init__(variational_strategy, input_dims, output_dims)
37 
38         self.num_quad_sites = num_quad_sites
39 
40         # Pass in previous_layer.quad_sites if you want to share quad_sites across layers.
41         if quad_sites is not None:
42             self.quad_sites = quad_sites
43         else:
44             self.quad_sites = torch.nn.Parameter(torch.randn(num_quad_sites, input_dims))
45 
46     def __call__(self, inputs, are_samples=False, expand_for_quadgrid=True, **kwargs):
47         if isinstance(inputs, MultitaskMultivariateNormal):
48             # inputs is definitely in the second layer, and mean is n x t
49             mus, sigmas = inputs.mean, inputs.variance.sqrt()
50 
51             if expand_for_quadgrid:
52                 xi_mus = mus.unsqueeze(0)  # 1 x n x t
53                 xi_sigmas = sigmas.unsqueeze(0)  # 1 x n x t
54             else:
55                 xi_mus = mus
56                 xi_sigmas = sigmas
57 
58             # unsqueeze sigmas to 1 x n x t, locations from [q] to Q^T x 1 x T.
59             # Broadcasted result will be Q^T x N x T
60             qg = self.quad_sites.view([self.num_quad_sites] + [1] * (xi_mus.dim() - 2) + [self.input_dims])
61             xi_sigmas = xi_sigmas * qg
62 
63             inputs = xi_mus + xi_sigmas  # q^t x n x t
64         if settings.debug.on():
65             if not torch.is_tensor(inputs):
66                 raise ValueError(
67                     "`inputs` should either be a MultitaskMultivariateNormal or a Tensor, got "
68                     f"{inputs.__class__.__Name__}"
69                 )
70 
71             if inputs.size(-1) != self.input_dims:
72                 raise RuntimeError(
73                     f"Input shape did not match self.input_dims. Got total feature dims [{inputs.size(-1)}],"
74                     f" expected [{self.input_dims}]"
75                 )
76 
77         # Repeat the input for all possible outputs
78         if self.output_dims is not None:
79             inputs = inputs.unsqueeze(-3)
80             inputs = inputs.expand(*inputs.shape[:-3], self.output_dims, *inputs.shape[-2:])
81         # Now run samples through the GP
82         output = ApproximateGP.__call__(self, inputs, **kwargs)
83 
84         if self.num_quad_sites > 0:
85             if self.output_dims is not None and not isinstance(output, MultitaskMultivariateNormal):
86                 mean = output.loc.transpose(-1, -2)
87                 covar = BlockDiagLazyTensor(output.lazy_covariance_matrix, block_dim=-3)
88                 output = MultitaskMultivariateNormal(mean, covar, interleaved=False)
89         else:
90             output = output.loc.transpose(-1, -2)  # this layer provides noiseless kernel interpolation
91 
92         return output
93 
94 
95 class DSPP(DeepGP):
96     """
97     A container module to build a DSPP
98     This module should contain :obj:`~gpytorch.models.deep_gps.DSPPLayer`
99     modules, and can also contain other modules as well.
100 
101     This Module contains an additional set of parameters, `raw_quad_weights`, that represent the mixture weights for
102     the output distribution.
103     """
104 
105     def __init__(self, num_quad_sites):
106         super().__init__()
107         self.num_quad_sites = num_quad_sites
108         self.register_parameter("raw_quad_weights", torch.nn.Parameter(torch.randn(self.num_quad_sites)))
109 
110     @property
111     def quad_weights(self):
112         qwd = self.raw_quad_weights
113         return qwd - qwd.logsumexp(dim=-1)
114
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/gpytorch/models/deep_gps/dspp.py b/gpytorch/models/deep_gps/dspp.py
--- a/gpytorch/models/deep_gps/dspp.py
+++ b/gpytorch/models/deep_gps/dspp.py
@@ -43,24 +43,18 @@
         else:
             self.quad_sites = torch.nn.Parameter(torch.randn(num_quad_sites, input_dims))
 
-    def __call__(self, inputs, are_samples=False, expand_for_quadgrid=True, **kwargs):
+    def __call__(self, inputs, **kwargs):
         if isinstance(inputs, MultitaskMultivariateNormal):
-            # inputs is definitely in the second layer, and mean is n x t
+            # This is for subsequent layers. We apply quadrature here
+            # Mean, stdv are q x ... x n x t
             mus, sigmas = inputs.mean, inputs.variance.sqrt()
+            qg = self.quad_sites.view([self.num_quad_sites] + [1] * (mus.dim() - 2) + [self.input_dims])
+            sigmas = sigmas * qg
+            inputs = mus + sigmas  # q^t x n x t
+            deterministic_inputs = False
+        else:
+            deterministic_inputs = True
 
-            if expand_for_quadgrid:
-                xi_mus = mus.unsqueeze(0)  # 1 x n x t
-                xi_sigmas = sigmas.unsqueeze(0)  # 1 x n x t
-            else:
-                xi_mus = mus
-                xi_sigmas = sigmas
-
-            # unsqueeze sigmas to 1 x n x t, locations from [q] to Q^T x 1 x T.
-            # Broadcasted result will be Q^T x N x T
-            qg = self.quad_sites.view([self.num_quad_sites] + [1] * (xi_mus.dim() - 2) + [self.input_dims])
-            xi_sigmas = xi_sigmas * qg
-
-            inputs = xi_mus + xi_sigmas  # q^t x n x t
         if settings.debug.on():
             if not torch.is_tensor(inputs):
                 raise ValueError(
@@ -78,9 +72,15 @@
         if self.output_dims is not None:
             inputs = inputs.unsqueeze(-3)
             inputs = inputs.expand(*inputs.shape[:-3], self.output_dims, *inputs.shape[-2:])
+
         # Now run samples through the GP
         output = ApproximateGP.__call__(self, inputs, **kwargs)
 
+        # If this is the first layer (deterministic inputs), expand the output
+        # This allows quadrature to be applied to future layers
+        if deterministic_inputs:
+            output = output.expand(torch.Size([self.num_quad_sites]) + output.batch_shape)
+
         if self.num_quad_sites > 0:
             if self.output_dims is not None and not isinstance(output, MultitaskMultivariateNormal):
                 mean = output.loc.transpose(-1, -2)

{"golden_diff": "diff --git a/gpytorch/models/deep_gps/dspp.py b/gpytorch/models/deep_gps/dspp.py\n--- a/gpytorch/models/deep_gps/dspp.py\n+++ b/gpytorch/models/deep_gps/dspp.py\n@@ -43,24 +43,18 @@\n         else:\n             self.quad_sites = torch.nn.Parameter(torch.randn(num_quad_sites, input_dims))\n \n-    def __call__(self, inputs, are_samples=False, expand_for_quadgrid=True, **kwargs):\n+    def __call__(self, inputs, **kwargs):\n         if isinstance(inputs, MultitaskMultivariateNormal):\n-            # inputs is definitely in the second layer, and mean is n x t\n+            # This is for subsequent layers. We apply quadrature here\n+            # Mean, stdv are q x ... x n x t\n             mus, sigmas = inputs.mean, inputs.variance.sqrt()\n+            qg = self.quad_sites.view([self.num_quad_sites] + [1] * (mus.dim() - 2) + [self.input_dims])\n+            sigmas = sigmas * qg\n+            inputs = mus + sigmas  # q^t x n x t\n+            deterministic_inputs = False\n+        else:\n+            deterministic_inputs = True\n \n-            if expand_for_quadgrid:\n-                xi_mus = mus.unsqueeze(0)  # 1 x n x t\n-                xi_sigmas = sigmas.unsqueeze(0)  # 1 x n x t\n-            else:\n-                xi_mus = mus\n-                xi_sigmas = sigmas\n-\n-            # unsqueeze sigmas to 1 x n x t, locations from [q] to Q^T x 1 x T.\n-            # Broadcasted result will be Q^T x N x T\n-            qg = self.quad_sites.view([self.num_quad_sites] + [1] * (xi_mus.dim() - 2) + [self.input_dims])\n-            xi_sigmas = xi_sigmas * qg\n-\n-            inputs = xi_mus + xi_sigmas  # q^t x n x t\n         if settings.debug.on():\n             if not torch.is_tensor(inputs):\n                 raise ValueError(\n@@ -78,9 +72,15 @@\n         if self.output_dims is not None:\n             inputs = inputs.unsqueeze(-3)\n             inputs = inputs.expand(*inputs.shape[:-3], self.output_dims, *inputs.shape[-2:])\n+\n         # Now run samples through the GP\n         output = ApproximateGP.__call__(self, inputs, **kwargs)\n \n+        # If this is the first layer (deterministic inputs), expand the output\n+        # This allows quadrature to be applied to future layers\n+        if deterministic_inputs:\n+            output = output.expand(torch.Size([self.num_quad_sites]) + output.batch_shape)\n+\n         if self.num_quad_sites > 0:\n             if self.output_dims is not None and not isinstance(output, MultitaskMultivariateNormal):\n                 mean = output.loc.transpose(-1, -2)\n", "issue": "[Docs] Deep Sigma Point Processes\nHi,\r\n\r\nThe DSPPs documentaiton makes it seem like theyare coded up very similarly to the standard Deep GPs, and this does seem true for[ two layers](https://docs.gpytorch.ai/en/stable/examples/05_Deep_Gaussian_Processes/Deep_Sigma_Point_Processes.html). \r\n\r\nHowever, if I try to add a third layer (a second hidden layer) and update the input and output dims and make sure the forward call is correct, the loss function returns a vector rather than a single value (which the optimsier obviously doesn't like).\r\n\r\nI'm guessing these losses are meant to be weighted according to the quadrature points. Is this true? Perhpas this could be more clear in the documentation (or alternatively, should this be done automatically in the backend?)\n", "before_files": [{"content": "import torch\n\nfrom gpytorch import settings\nfrom gpytorch.distributions import MultitaskMultivariateNormal\nfrom gpytorch.lazy import BlockDiagLazyTensor\n\nfrom ..approximate_gp import ApproximateGP\nfrom .deep_gp import DeepGP, DeepGPLayer\n\n\nclass DSPPLayer(DeepGPLayer):\n    \"\"\"\n    Represents a layer in a DSPP where inference is performed using the techniques defined in Jankowiak et al., 2020.\n\n    From an end user perspective, the functionality and usage of a DSPPLayer is essentially identical to that of a\n    DeepGPLayer. It is therefore recommended that you review the documentation for DeepGPLayer.\n\n    :param ~gpytorch.variational.VariationalStrategy variational_strategy: Strategy for\n        changing q(u) -> q(f) (see other VI docs)\n    :param int input_dims: Dimensionality of input data expected by each GP\n    :param int output_dims: (default None) Number of GPs in this layer, equivalent to\n        output dimensionality. If set to `None`, then the output dimension will be squashed.\n    :param int num_quad_sites: Number of quadrature sites to use. Also the number of Gaussians in the mixture output\n        by this layer.\n\n    Again, refer to the documentation for DeepGPLayer or our example notebooks for full details on what calling a\n    DSPPLayer module does. The high level overview is that if a tensor `x` is `n x d` then\n\n        >>> hidden_gp2(hidden_gp1(x))\n\n    will return a `num_quad_sites` by `output_dims` set of Gaussians, where for each output dim the first batch dim\n    represents a weighted mixture of `num_quad_sites` Gaussians with weights given by DSPP.quad_weights (see DSPP below)\n    \"\"\"\n\n    def __init__(self, variational_strategy, input_dims, output_dims, num_quad_sites=3, quad_sites=None):\n        super().__init__(variational_strategy, input_dims, output_dims)\n\n        self.num_quad_sites = num_quad_sites\n\n        # Pass in previous_layer.quad_sites if you want to share quad_sites across layers.\n        if quad_sites is not None:\n            self.quad_sites = quad_sites\n        else:\n            self.quad_sites = torch.nn.Parameter(torch.randn(num_quad_sites, input_dims))\n\n    def __call__(self, inputs, are_samples=False, expand_for_quadgrid=True, **kwargs):\n        if isinstance(inputs, MultitaskMultivariateNormal):\n            # inputs is definitely in the second layer, and mean is n x t\n            mus, sigmas = inputs.mean, inputs.variance.sqrt()\n\n            if expand_for_quadgrid:\n                xi_mus = mus.unsqueeze(0)  # 1 x n x t\n                xi_sigmas = sigmas.unsqueeze(0)  # 1 x n x t\n            else:\n                xi_mus = mus\n                xi_sigmas = sigmas\n\n            # unsqueeze sigmas to 1 x n x t, locations from [q] to Q^T x 1 x T.\n            # Broadcasted result will be Q^T x N x T\n            qg = self.quad_sites.view([self.num_quad_sites] + [1] * (xi_mus.dim() - 2) + [self.input_dims])\n            xi_sigmas = xi_sigmas * qg\n\n            inputs = xi_mus + xi_sigmas  # q^t x n x t\n        if settings.debug.on():\n            if not torch.is_tensor(inputs):\n                raise ValueError(\n                    \"`inputs` should either be a MultitaskMultivariateNormal or a Tensor, got \"\n                    f\"{inputs.__class__.__Name__}\"\n                )\n\n            if inputs.size(-1) != self.input_dims:\n                raise RuntimeError(\n                    f\"Input shape did not match self.input_dims. Got total feature dims [{inputs.size(-1)}],\"\n                    f\" expected [{self.input_dims}]\"\n                )\n\n        # Repeat the input for all possible outputs\n        if self.output_dims is not None:\n            inputs = inputs.unsqueeze(-3)\n            inputs = inputs.expand(*inputs.shape[:-3], self.output_dims, *inputs.shape[-2:])\n        # Now run samples through the GP\n        output = ApproximateGP.__call__(self, inputs, **kwargs)\n\n        if self.num_quad_sites > 0:\n            if self.output_dims is not None and not isinstance(output, MultitaskMultivariateNormal):\n                mean = output.loc.transpose(-1, -2)\n                covar = BlockDiagLazyTensor(output.lazy_covariance_matrix, block_dim=-3)\n                output = MultitaskMultivariateNormal(mean, covar, interleaved=False)\n        else:\n            output = output.loc.transpose(-1, -2)  # this layer provides noiseless kernel interpolation\n\n        return output\n\n\nclass DSPP(DeepGP):\n    \"\"\"\n    A container module to build a DSPP\n    This module should contain :obj:`~gpytorch.models.deep_gps.DSPPLayer`\n    modules, and can also contain other modules as well.\n\n    This Module contains an additional set of parameters, `raw_quad_weights`, that represent the mixture weights for\n    the output distribution.\n    \"\"\"\n\n    def __init__(self, num_quad_sites):\n        super().__init__()\n        self.num_quad_sites = num_quad_sites\n        self.register_parameter(\"raw_quad_weights\", torch.nn.Parameter(torch.randn(self.num_quad_sites)))\n\n    @property\n    def quad_weights(self):\n        qwd = self.raw_quad_weights\n        return qwd - qwd.logsumexp(dim=-1)\n", "path": "gpytorch/models/deep_gps/dspp.py"}], "after_files": [{"content": "import torch\n\nfrom gpytorch import settings\nfrom gpytorch.distributions import MultitaskMultivariateNormal\nfrom gpytorch.lazy import BlockDiagLazyTensor\n\nfrom ..approximate_gp import ApproximateGP\nfrom .deep_gp import DeepGP, DeepGPLayer\n\n\nclass DSPPLayer(DeepGPLayer):\n    \"\"\"\n    Represents a layer in a DSPP where inference is performed using the techniques defined in Jankowiak et al., 2020.\n\n    From an end user perspective, the functionality and usage of a DSPPLayer is essentially identical to that of a\n    DeepGPLayer. It is therefore recommended that you review the documentation for DeepGPLayer.\n\n    :param ~gpytorch.variational.VariationalStrategy variational_strategy: Strategy for\n        changing q(u) -> q(f) (see other VI docs)\n    :param int input_dims: Dimensionality of input data expected by each GP\n    :param int output_dims: (default None) Number of GPs in this layer, equivalent to\n        output dimensionality. If set to `None`, then the output dimension will be squashed.\n    :param int num_quad_sites: Number of quadrature sites to use. Also the number of Gaussians in the mixture output\n        by this layer.\n\n    Again, refer to the documentation for DeepGPLayer or our example notebooks for full details on what calling a\n    DSPPLayer module does. The high level overview is that if a tensor `x` is `n x d` then\n\n        >>> hidden_gp2(hidden_gp1(x))\n\n    will return a `num_quad_sites` by `output_dims` set of Gaussians, where for each output dim the first batch dim\n    represents a weighted mixture of `num_quad_sites` Gaussians with weights given by DSPP.quad_weights (see DSPP below)\n    \"\"\"\n\n    def __init__(self, variational_strategy, input_dims, output_dims, num_quad_sites=3, quad_sites=None):\n        super().__init__(variational_strategy, input_dims, output_dims)\n\n        self.num_quad_sites = num_quad_sites\n\n        # Pass in previous_layer.quad_sites if you want to share quad_sites across layers.\n        if quad_sites is not None:\n            self.quad_sites = quad_sites\n        else:\n            self.quad_sites = torch.nn.Parameter(torch.randn(num_quad_sites, input_dims))\n\n    def __call__(self, inputs, **kwargs):\n        if isinstance(inputs, MultitaskMultivariateNormal):\n            # This is for subsequent layers. We apply quadrature here\n            # Mean, stdv are q x ... x n x t\n            mus, sigmas = inputs.mean, inputs.variance.sqrt()\n            qg = self.quad_sites.view([self.num_quad_sites] + [1] * (mus.dim() - 2) + [self.input_dims])\n            sigmas = sigmas * qg\n            inputs = mus + sigmas  # q^t x n x t\n            deterministic_inputs = False\n        else:\n            deterministic_inputs = True\n\n        if settings.debug.on():\n            if not torch.is_tensor(inputs):\n                raise ValueError(\n                    \"`inputs` should either be a MultitaskMultivariateNormal or a Tensor, got \"\n                    f\"{inputs.__class__.__Name__}\"\n                )\n\n            if inputs.size(-1) != self.input_dims:\n                raise RuntimeError(\n                    f\"Input shape did not match self.input_dims. Got total feature dims [{inputs.size(-1)}],\"\n                    f\" expected [{self.input_dims}]\"\n                )\n\n        # Repeat the input for all possible outputs\n        if self.output_dims is not None:\n            inputs = inputs.unsqueeze(-3)\n            inputs = inputs.expand(*inputs.shape[:-3], self.output_dims, *inputs.shape[-2:])\n\n        # Now run samples through the GP\n        output = ApproximateGP.__call__(self, inputs, **kwargs)\n\n        # If this is the first layer (deterministic inputs), expand the output\n        # This allows quadrature to be applied to future layers\n        if deterministic_inputs:\n            output = output.expand(torch.Size([self.num_quad_sites]) + output.batch_shape)\n\n        if self.num_quad_sites > 0:\n            if self.output_dims is not None and not isinstance(output, MultitaskMultivariateNormal):\n                mean = output.loc.transpose(-1, -2)\n                covar = BlockDiagLazyTensor(output.lazy_covariance_matrix, block_dim=-3)\n                output = MultitaskMultivariateNormal(mean, covar, interleaved=False)\n        else:\n            output = output.loc.transpose(-1, -2)  # this layer provides noiseless kernel interpolation\n\n        return output\n\n\nclass DSPP(DeepGP):\n    \"\"\"\n    A container module to build a DSPP\n    This module should contain :obj:`~gpytorch.models.deep_gps.DSPPLayer`\n    modules, and can also contain other modules as well.\n\n    This Module contains an additional set of parameters, `raw_quad_weights`, that represent the mixture weights for\n    the output distribution.\n    \"\"\"\n\n    def __init__(self, num_quad_sites):\n        super().__init__()\n        self.num_quad_sites = num_quad_sites\n        self.register_parameter(\"raw_quad_weights\", torch.nn.Parameter(torch.randn(self.num_quad_sites)))\n\n    @property\n    def quad_weights(self):\n        qwd = self.raw_quad_weights\n        return qwd - qwd.logsumexp(dim=-1)\n", "path": "gpytorch/models/deep_gps/dspp.py"}]}