Return the predictive mean and variance at a series of new points X1new, X2new
        Only returns the diagonal of the predictive variance, for now.

        :param X1new: The points at which to make a prediction
        :type X1new: np.ndarray, Nnew x self.input_dim1
        :param X2new: The points at which to make a prediction
        :type X2new: np.ndarray, Nnew x self.input_dim2

def predict(self, X1new, X2new):
        """
        Return the predictive mean and variance at a series of new points X1new, X2new
        Only returns the diagonal of the predictive variance, for now.

        :param X1new: The points at which to make a prediction
        :type X1new: np.ndarray, Nnew x self.input_dim1
        :param X2new: The points at which to make a prediction
        :type X2new: np.ndarray, Nnew x self.input_dim2

        """
        k1xf = self.kern1.K(X1new, self.X1)
        k2xf = self.kern2.K(X2new, self.X2)
        A = k1xf.dot(self.U1)
        B = k2xf.dot(self.U2)
        mu = A.dot(self.Ytilde.reshape(self.num_data1, self.num_data2, order='F')).dot(B.T).flatten(order='F')
        k1xx = self.kern1.Kdiag(X1new)
        k2xx = self.kern2.Kdiag(X2new)
        BA = np.kron(B, A)
        var = np.kron(k2xx, k1xx) - np.sum(BA**2*self.Wi, 1) + self.likelihood.variance

        return mu[:, None], var[:, None]