code: def _fD1(self, neg_pairs, A):
      
      dim = neg_pairs.shape[2]
      diff = neg_pairs[:, 0, :] - neg_pairs[:, 1, :]
      # outer products of all rows in `diff`
      M = np.einsum('ij,ik->ijk', diff, diff)
      # faster version of: dist = np.sqrt(np.sum(M * A[None,:,:], axis=(1,2)))
      dist = np.sqrt(np.einsum('ijk,jk', M, A))
      # faster version of: sum_deri = np.sum(M / (2 * (dist[:,None,None] + 1e-6)), axis=0)
      sum_deri = np.einsum('ijk,i->jk', M, 0.5 / (dist + 1e-6))
      sum_dist = dist.sum()
      return sum_deri / (sum_dist + 1e-6)
docs: The gradient of the dissimilarity constraint function w.r.t. A.
  
      For example, let distance by L1 norm:
      f = f(\sum_{ij \in D} \sqrt{(x_i-x_j)A(x_i-x_j)'})
      df/dA_{kl} = f'* d(\sum_{ij \in D} \sqrt{(x_i-x_j)^k*(x_i-x_j)^l})/dA_{kl}
  
      Note that d_ij*A*d_ij' = tr(d_ij*A*d_ij') = tr(d_ij'*d_ij*A)
      so, d(d_ij*A*d_ij')/dA = d_ij'*d_ij
          df/dA = f'(\sum_{ij \in D} \sqrt{tr(d_ij'*d_ij*A)})
                  * 0.5*(\sum_{ij \in D} (1/sqrt{tr(d_ij'*d_ij*A)})*(d_ij'*d_ij))