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
Find Sopt given X, Y