code: def _expectation(p, kern1, feat1, kern2, feat2, nghp=None):
      r
      crossexps = []
  
      if kern1 == kern2 and feat1 == feat2:  # avoid duplicate computation by using transposes
          for i, k1 in enumerate(kern1.kernels):
              crossexps.append(expectation(p, (k1, feat1), (k1, feat1), nghp=nghp))
  
              for k2 in kern1.kernels[:i]:
                  eKK = expectation(p, (k1, feat1), (k2, feat2), nghp=nghp)
                  eKK += tf.matrix_transpose(eKK)
                  crossexps.append(eKK)
      else:
          for k1, k2 in it.product(kern1.kernels, kern2.kernels):
              crossexps.append(expectation(p, (k1, feat1), (k2, feat2), nghp=nghp))
  
      return functools.reduce(tf.add, crossexps)
docs: r"""
      Compute the expectation:
      expectation[n] = <(\Sum_i K1_i_{Z1, x_n}) (\Sum_j K2_j_{x_n, Z2})>_p(x_n)
          - \Sum_i K1_i_{.,.}, \Sum_j K2_j_{.,.} :: Sum kernels
  
      :return: NxM1xM2