def predict(self, X1new, X2new):
k1xf = self.kern1.K(X1new, self.X1)<EOL>k2xf = self.kern2.K(X2new, self.X2)<EOL>A = k1xf.dot(self.U1)<EOL>B = k2xf.dot(self.U2)<EOL>mu = A.dot(self.Ytilde.reshape(self.num_data1, self.num_data2, order='<STR_LIT:F>')).dot(B.T).flatten(order='<STR_LIT:F>')<EOL>k1xx = self.kern1.Kdiag(X1new)<EOL>k2xx = self.kern2.Kdiag(X2new)<EOL>BA = np.kron(B, A)<EOL>var = np.kron(k2xx, k1xx) - np.sum(BA**<NUM_LIT:2>*self.Wi, <NUM_LIT:1>) + self.likelihood.variance<EOL>return mu[:, None], var[:, None]<EOL>
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