code: def _fS1(self, pos_pairs, A):
      
      dim = pos_pairs.shape[2]
      diff = pos_pairs[:, 0, :] - pos_pairs[:, 1, :]
      return np.einsum('ij,ik->jk', diff, diff)
docs: The gradient of the similarity constraint function w.r.t. A.
  
      f = \sum_{ij}(x_i-x_j)A(x_i-x_j)' = \sum_{ij}d_ij*A*d_ij'
      df/dA = d(d_ij*A*d_ij')/dA
  
      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