def METHOD_NAME(
    X: tf.Tensor, Z: tf.Tensor, kern: Kernel, q_mu: tf.Tensor, q_sqrt: tf.Tensor
) -> MeanAndVariance:
    Kmm = kern(Z, Z) + np.eye(M) * gpflow.config.default_jitter()
    Kmn = kern(Z, X)
    Knn = kern(X, X)
    Kmm, Kmn, Knn = [k.numpy() for k in [Kmm, Kmn, Knn]]
    Knm: AnyNDArray = Kmn.T
    Kmm, Kmn, Knm, Knn = [np.tile(k[None, :, :], [Dy, 1, 1]) for k in [Kmm, Kmn, Knm, Knn]]
    S = q_sqrt @ np.transpose(q_sqrt, [0, 2, 1])
    Kmm_inv = np.linalg.inv(Kmm)
    mean = np.einsum("dmn,dmM,Md->nd", Kmn, Kmm_inv, q_mu)
    cov = Knn + Knm @ Kmm_inv @ (S - Kmm) @ Kmm_inv @ Kmn
    return mean, cov