docs: def _expectation(p, kern, feat, mean, none, nghp=None):
      """
      Compute the expectation:
      expectation[n] = <K_{Z, x_n} x_{n+1}^T>_p(x_{n:n+1})
          - K_{.,.} :: Linear kernel
          - p       :: MarkovGaussian distribution (p.cov 2x(N+1)xDxD)
  
      :return: NxMxD
      """
      Xmu, Xcov = p.mu, p.cov
  
      with tf.control_dependencies([tf.assert_equal(
              tf.shape(Xmu)[1], tf.constant(kern.input_dim, settings.int_type),
              message="Currently cannot handle slicing in exKxz.")]):
          Xmu = tf.identity(Xmu)
  
      with params_as_tensors_for(kern, feat):
          N = tf.shape(Xmu)[0] - 1
          var_Z = kern.variance * feat.Z  # MxD
          tiled_Z = tf.tile(tf.expand_dims(var_Z, 0), (N, 1, 1))  # NxMxD
          eXX = Xcov[1, :-1] + (Xmu[:-1][..., None] * Xmu[1:][:, None, :])  # NxDxD
          return tf.matmul(tiled_Z, eXX)