def power_method(U_hat, U, tol=1e-4):
    n, _ = U.shape
    x = np.random.normal(size=(n,))
    s = np.linalg.norm(x)
    x = x / s
    while True:
        x = np.dot(U_hat, np.dot(U_hat.T, x)) - np.dot(U, np.dot(U.T, x))
        s_new = np.linalg.norm(x)
        if np.abs(s - s_new) < tol:
            break
        s = s_new
        x = x / s
    return s