docs: def independent_interdomain_conditional(Kmn, Kmm, Knn, f, *, full_cov=False, full_output_cov=False,
                                          q_sqrt=None, white=False):
      """
      The inducing outputs live in the g-space (R^L).
      Interdomain conditional calculation.
  
      :param Kmn: M x L x N x P
      :param Kmm: L x M x M
      :param Knn: N x P  or  N x N  or  P x N x N  or  N x P x N x P
      :param f: data matrix, M x L
      :param q_sqrt: L x M x M  or  M x L
      :param full_cov: calculate covariance between inputs
      :param full_output_cov: calculate covariance between outputs
      :param white: use whitened representation
      :return:
          - mean: N x P
          - variance: N x P, N x P x P, P x N x N, N x P x N x P
      """
      logger.debug("independent_interdomain_conditional")
      M, L, N, P = [tf.shape(Kmn)[i] for i in range(Kmn.shape.ndims)]
  
      Lm = tf.cholesky(Kmm)  # L x M x M
  
      # Compute the projection matrix A
      Kmn = tf.reshape(tf.transpose(Kmn, (1, 0, 2, 3)), (L, M, N * P))
      A = tf.matrix_triangular_solve(Lm, Kmn, lower=True)  # L x M x M  *  L x M x NP  ->  L x M x NP
      Ar = tf.reshape(A, (L, M, N, P))
  
      # compute the covariance due to the conditioning
      if full_cov and full_output_cov:
          fvar = Knn - tf.tensordot(Ar, Ar, [[0, 1], [0, 1]])  # N x P x N x P
      elif full_cov and not full_output_cov:
          At = tf.reshape(tf.transpose(Ar), (P, N, M * L))  # P x N x ML
          fvar = Knn - tf.matmul(At, At, transpose_b=True)  # P x N x N
      elif not full_cov and full_output_cov:
          At = tf.reshape(tf.transpose(Ar, [2, 3, 1, 0]), (N, P, M * L))  # N x P x ML
          fvar = Knn - tf.matmul(At, At, transpose_b=True)  # N x P x P
      elif not full_cov and not full_output_cov:
          fvar = Knn - tf.reshape(tf.reduce_sum(tf.square(A), [0, 1]), (N, P))  # Knn: N x P
  
      # another backsubstitution in the unwhitened case
      if not white:
          A = tf.matrix_triangular_solve(Lm, Ar)  # L x M x M  *  L x M x NP  ->  L x M x NP
          Ar = tf.reshape(A, (L, M, N, P))
  
      fmean = tf.tensordot(Ar, f, [[1, 0], [0, 1]])  # N x P
  
      if q_sqrt is not None:
          if q_sqrt.shape.ndims == 3:
              Lf = tf.matrix_band_part(q_sqrt, -1, 0)  # L x M x M
              LTA = tf.matmul(Lf, A, transpose_a=True)  # L x M x M  *  L x M x NP  ->  L x M x NP
          else:  # q_sqrt M x L
              LTA = (A * tf.transpose(q_sqrt)[..., None])  # L x M x NP
  
          if full_cov and full_output_cov:
              LTAr = tf.reshape(LTA, (L * M, N * P))
              fvar = fvar + tf.reshape(tf.matmul(LTAr, LTAr, transpose_a=True), (N, P, N, P))
          elif full_cov and not full_output_cov:
              LTAr = tf.transpose(tf.reshape(LTA, (L * M, N, P)), [2, 0, 1])  # P x LM x N
              fvar = fvar + tf.matmul(LTAr, LTAr, transpose_a=True)  # P x N x N
          elif not full_cov and full_output_cov:
              LTAr = tf.transpose(tf.reshape(LTA, (L * M, N, P)), [1, 0, 2])  # N x LM x P
              fvar = fvar + tf.matmul(LTAr, LTAr, transpose_a=True)  # N x P x P
          elif not full_cov and not full_output_cov:
              fvar = fvar + tf.reshape(tf.reduce_sum(tf.square(LTA), (0, 1)), (N, P))
      return fmean, fvar