docs: def getoptS(X, Y, M_E, E):
      ''' Find Sopt given X, Y
      '''
      n, r = X.shape
      C = np.dot(np.dot(X.T, M_E), Y)
      C = C.flatten()
  
      A = np.zeros((r * r, r * r))
      for i in range(r):
          for j in range(r):
              ind = j * r + i
              temp = np.dot(
                  np.dot(X.T, np.dot(X[:, i, None], Y[:, j, None].T) * E), Y)
              A[:, ind] = temp.flatten()
  
      S = np.linalg.solve(A, C)
  
      return np.reshape(S, (r, r)).T