[Bug] Deep Kernel Learning -  face orientation extraction
# 🐛 Bug

Hi,

I am trying to implement the olivetti face orientation classification problem from https://arxiv.org/abs/1511.02222 with gpytorch. Unfortunetly I could not manage to solve the following error.

## To reproduce

```python
import os
import tqdm
import torch
import gpytorch
import numpy as np
from matplotlib import pyplot as plt


# Define a plotting function
def ax_plot(f, ax, y_labels, title):
    im = ax.imshow(y_labels)
    ax.set_title(title)
    f.colorbar(im)
    

# Load Olivetti faces database
data = np.load('olivetti.npz')

train_x = torch.Tensor(data['train_x'])
train_x = train_x.reshape((train_x.shape[0], 1, train_x.shape[1], train_x.shape[2]))
train_y = torch.Tensor(data['train_y'])


test_x = torch.Tensor(data['test_x'])
test_x = test_x.reshape((test_x.shape[0], 1, test_x.shape[1], test_x.shape[2]))
test_y = torch.Tensor(data['test_y'])


# Create dataset
dataset    = torch.utils.data.TensorDataset(train_x, train_y)
dataloader = torch.utils.data.DataLoader(dataset) 

if torch.cuda.is_available():
    train_x, train_y = train_x.cuda(), train_y.cuda()
    
    
# Defining the DKL Feature Extractor
data_dim = train_x.size(-1)

class LargeFeatureExtractor(torch.nn.Sequential):           
    def __init__(self):                                      
        super(LargeFeatureExtractor, self).__init__()        
        self.add_module('conv1', torch.nn.Conv2d(kernel_size=(5,5), stride=1, in_channels=1, out_channels=20))
        self.add_module('pool1', torch.nn.MaxPool2d(kernel_size=(2, 2), stride=2))                  
        self.add_module('conv2', torch.nn.Conv2d(kernel_size=(5,5), stride=1, in_channels=20, out_channels=50))     
        self.add_module('pool2', torch.nn.MaxPool2d(kernel_size=(2, 2), stride=2))
        self.add_module('flatt', torch.nn.Flatten())              
        self.add_module('full3', torch.nn.Linear(800, 1000))
        self.add_module('relu3', torch.nn.ReLU())
        self.add_module('full4', torch.nn.Linear(1000, 500))                  
        self.add_module('relu4', torch.nn.ReLU())
        self.add_module('full5', torch.nn.Linear(500, 50))       
        self.add_module('full6', torch.nn.Linear(50, 2))
        
feature_extractor = LargeFeatureExtractor()


# Defining the DKL-GP Model
class GPRegressionModel(gpytorch.models.ExactGP):
        def __init__(self, train_x, train_y, likelihood):
            super(GPRegressionModel, self).__init__(train_x, train_y, likelihood)
            self.mean_module = gpytorch.means.ConstantMean()
            self.covar_module = gpytorch.kernels.GridInterpolationKernel(
                gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel(ard_num_dims=1)),
                num_dims=2, grid_size=100
            )
            self.feature_extractor = feature_extractor

        def forward(self, x):
            # We're first putting our data through a deep net (feature extractor)
            # We're also scaling the features so that they're nice values
            projected_x = self.feature_extractor(x)
            projected_x = projected_x - projected_x.min(0)[0]
            projected_x = 2 * (projected_x / projected_x.max(0)[0]) - 1
        
            mean_x = self.mean_module(projected_x)
            covar_x = self.covar_module(projected_x)
            return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

        

likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = GPRegressionModel(train_x, train_y, likelihood)

if torch.cuda.is_available():
    model = model.cuda()
    likelihood = likelihood.cuda()
    
# Training the model
training_iterations = 1

# Find optimal model hyperparameters
model.train()
likelihood.train()

# Use the adam optimizer
optimizer = torch.optim.Adam([
    {'params': model.feature_extractor.parameters()},
    {'params': model.covar_module.parameters()},
    {'params': model.mean_module.parameters()},
    {'params': model.likelihood.parameters()},
], lr=0.01)

# "Loss" for GPs - the marginal log likelihood
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

def train():
    iterator = tqdm.tqdm_notebook(range(training_iterations))
    for i in iterator:
        # Zero backprop gradients
        optimizer.zero_grad()
        # Get output from model
        output = model(train_x)
        # Calc loss and backprop derivatives
        loss = -mll(output, train_y)
        loss.backward()
        iterator.set_postfix(loss=loss.item())
        optimizer.step()
        
# train model
train()


# Making predictions
model.eval()
likelihood.eval()
with torch.no_grad(), gpytorch.settings.use_toeplitz(False), gpytorch.settings.fast_pred_var():
    preds = model(test_x)

print('Test MAE: {}'.format(torch.mean(torch.abs(preds.mean - test_y))))
```

** Stack trace/error message **
```
Traceback (most recent call last):

  File "<ipython-input-1-88b9df752e72>", line 1, in <module>
    runfile('/home/user/Desktop/dl-dev python/gp_examples/test.py', wdir='/home/user/Desktop/dl-dev python/gp_examples')

  File "/home/user/anaconda3/envs/dkl/lib/python3.6/site-packages/spyder_kernels/customize/spydercustomize.py", line 827, in runfile
    execfile(filename, namespace)

  File "/home/user/anaconda3/envs/dkl/lib/python3.6/site-packages/spyder_kernels/customize/spydercustomize.py", line 110, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "/home/user/Desktop/dl-dev python/gp_examples/test.py", line 136, in <module>
    preds = model(test_x)

  File "/home/user/anaconda3/envs/dkl/lib/python3.6/site-packages/gpytorch/models/exact_gp.py", line 294, in __call__
    likelihood=self.likelihood,

  File "/home/user/anaconda3/envs/dkl/lib/python3.6/site-packages/gpytorch/models/exact_prediction_strategies.py", line 36, in prediction_strategy
    return cls(train_inputs, train_prior_dist, train_labels, likelihood)

  File "/home/user/anaconda3/envs/dkl/lib/python3.6/site-packages/gpytorch/kernels/grid_interpolation_kernel.py", line 205, in prediction_strategy
    return InterpolatedPredictionStrategy(train_inputs, train_prior_dist, train_labels, likelihood)

  File "/home/user/anaconda3/envs/dkl/lib/python3.6/site-packages/gpytorch/models/exact_prediction_strategies.py", line 380, in __init__
    super().__init__(train_inputs, train_prior_dist, train_labels, likelihood)

  File "/home/user/anaconda3/envs/dkl/lib/python3.6/site-packages/gpytorch/models/exact_prediction_strategies.py", line 50, in __init__
    mvn = self.likelihood(train_prior_dist, train_inputs)

  File "/home/user/anaconda3/envs/dkl/lib/python3.6/site-packages/gpytorch/likelihoods/likelihood.py", line 313, in __call__
    return self.marginal(input, *args, **kwargs)

  File "/home/user/anaconda3/envs/dkl/lib/python3.6/site-packages/gpytorch/likelihoods/gaussian_likelihood.py", line 76, in marginal
    full_covar = covar + noise_covar

  File "/home/user/anaconda3/envs/dkl/lib/python3.6/site-packages/gpytorch/lazy/lazy_tensor.py", line 1626, in __add__
    return AddedDiagLazyTensor(self, other)

  File "/home/user/anaconda3/envs/dkl/lib/python3.6/site-packages/gpytorch/lazy/added_diag_lazy_tensor.py", line 27, in __init__
    broadcasting._mul_broadcast_shape(lazy_tensors[0].shape, lazy_tensors[1].shape)

  File "/home/user/anaconda3/envs/dkl/lib/python3.6/site-packages/gpytorch/utils/broadcasting.py", line 20, in _mul_broadcast_shape
    raise RuntimeError("Shapes are not broadcastable for mul operation")

RuntimeError: Shapes are not broadcastable for mul operatio
```

## Expected Behavior

I seems that the train mode is working well, but the eval mode fails.

## System information
- GPyTorch Version 1.0.1
- PyTorch Version 1.4.0
- Computer OS Ubuntu 18.04.4 LTS

## Additional context
I may be missing something in the feature extractor...
[dkl.zip](https://github.com/cornellius-gp/gpytorch/files/4380811/dkl.zip)
diff --git a/gpytorch/likelihoods/gaussian_likelihood.py b/gpytorch/likelihoods/gaussian_likelihood.py
index 6080c887..affa636b 100644
--- a/gpytorch/likelihoods/gaussian_likelihood.py
+++ b/gpytorch/likelihoods/gaussian_likelihood.py
@@ -30,13 +30,7 @@ class _GaussianLikelihoodBase(Likelihood):
         self.noise_covar = noise_covar
 
     def _shaped_noise_covar(self, base_shape: torch.Size, *params: Any, **kwargs: Any):
-        if len(params) > 0:
-            # we can infer the shape from the params
-            shape = None
-        else:
-            # here shape[:-1] is the batch shape requested, and shape[-1] is `n`, the number of points
-            shape = base_shape
-        return self.noise_covar(*params, shape=shape, **kwargs)
+        return self.noise_covar(*params, shape=base_shape, **kwargs)
 
     def expected_log_prob(self, target: Tensor, input: MultivariateNormal, *params: Any, **kwargs: Any) -> Tensor:
         mean, variance = input.mean, input.variance
['gpytorch/likelihoods/gaussian_likelihood.py']