Solution 1:

Here's a bit of relevant literature:

  • Deep Sets (NIPS 2017) classifies all linear functions $\mathbb R^{n\times k} \to \mathbb R^{n\times l}$ that are permutation invariant / equivariant across the first axis.

  • On Universal Equivariant Set Networks (ICLR 2020) deals with the case of finding all homogeneous polynomial functions $\mathbb R^{n\times k} \to \mathbb R^{n\times l}$ that are permutation equivariant across the first axis.

  • Invariant and Equivariant Graph Networks (ICLR 2019) deals with the case of finding linear functions $\mathbb R^{n^k} \to \mathbb R^{n^k}$ that are permutation invariant / equivariant across all $k$ axes

Especially, the 3rd paper gives as an example the case of linear function $\mathbb R^{n^2}\to\mathbb R^{n^2}$ that are permutation equivariant across each axis, i.e. $f(P^T X P)=P^T f(X) P$, which is precisely your problem. They show that the space of such linear functions is $15$-dimensional, independent of $n$ (!).

One can combine papers 2 and 3 to find all homogeneous polynomial functions $\mathbb R^{n^k}\to\mathbb R^{n^k}$ that are permutation equivariant across all $k$ axes.