Question
Show that there exists a set of unitary matrices $\{U_i\}$, and probability $\{p_i\}$, such that for any $n \times n$ matrix $A$ \begin{equation} \tag{1} \sum_{i} p_i U_i A U^{\dagger}_i = \text{tr}(A) \frac{I}{n} \end{equation}
Attempts
For $n=2$, it is easy to show \begin{equation} \frac{1}{4} ( \sigma^x A \sigma^x + \sigma^y A \sigma^y + \sigma^z A \sigma^z + I A I ) = \text{tr}(A) I / 2 \end{equation} where $\sigma^{x,y,z}$ are Pauli sigma matrices. The idea comes from kraus operator sum representation.
We can then generalize to dimension $n = 2^m$, where $U_i$ can be taken as the tensor products of these basis, but not arbitrary dimension.
In indices, Eq.(1) is equivalent to \begin{equation} \sum_i p_i (U_i)_{ab} (U_i^*)_{dc} = \delta_{bc} \delta_{ad} / n \end{equation} This looks like the identity from the finite dimensional irreducible unitary representation of finite group, see Peter-Weyl theorem. But again this only works when group $G$ has irreducible representation at dimension $n$, and all $p_i$ in this case are equal.
I feel that "right proof" should not utilize these additional structures.