[Take 4.2, back to requiring only the PRP to be secret]
Improving on the line of thought in that other answerthat other answer, we will craft two efficiently computable functions $H_1$ and $H_2$ each accepting two arbitrary strings as input, with output a fixed-size string, say 512 bit; and (I conjecture) indistinguishable from a random function (under the assumption that some parameters are secret) except for this superset of the desired property: $$\forall(i,j)\in\{1,2\}^2,\forall(s_1,s_2,x), H_i(s_1,H_j(s_2,x))=H_j(s_2,H_i(s_1,x))$$
We'll use as building blocks three 512-bit PRFs $P_0()$, $P_1()$ and $P_2()$ accepting an arbitrary bit string as input, and one 512-bit PRP $E()$ with $D()$ its reverse function, such that $D(E(a))=a$.
We can build the PRFs as $P(m)=\mathtt{HMAC}_\mathtt{SHA512}(k,m)$ and three arbitrary distinct $k$. And by this famous result we can build the PRP as the four-rounds Feistel cipher with round functions $\mathtt{HMAC}_\mathtt{SHA256}(k,m)$ and four arbitrary distinct $k$.
For $i\in\{1,2\}$, define $H_i(s,x)$ as $$H_i(s,x)=E\big(\text{ }D(P_i(s))+D(x)\bmod{2^{512}}\text{ }\big)\text{ when }x\text{ is a 512-bit string,}$$ $$H_i(s,x)=E\big(\text{ }D(P_i(s))+P_0(x)\bmod{2^{512}}\text{ }\big)\text{ otherwise.}$$Notice that in the above, $P_0(x)$ is now added without going through $D$; this avoid that knowledge of $P_0$ would allows creating collisions of the form $H_i(s,P_0(x))=H_i(s,x)$
The desired property $H_i(s_1,H_j(s_2,x))=H_j(s_2,H_i(s_1,x))$ follows from commutativity and associativity of addition in $\mathbb Z_{2^{512}}$, and the $H_i$ functions appears much like random oracles, except for consequences of that property, to an adversary ignoring the PRP (but possibly knowing the PRFs).
By querying the functions/oracles $H_i$, it seems impossible to internally add 0, or otherwise create collisions, other than by chance. Computing $H_1(s,H_1(s,\cdots H_1(s,H_1(s,x))\cdots))$, where $x$ is a 512-bit string and there are $2^n$ iterations, has odds only $2^{n-511-k}$ to return to $x$; if we had used XOR instead of modular addition, we'd have $H_1(s,H_1(s,x))=x$.
Note: to an adversary knowing the PRP and PRF, the problem of finding collisions on $H_1$ unrelated to the property is tractable: it is a knapsack problem in $\mathbb Z_{2^{512}}$ with unlimited supply of random values to choose from.
Open problems:
- Can this be proven secure (e.g. collision-resistant except per application of the required property); or broken?
- Can we make a scheme secure yet fully public?
- We do not need a full group operation (on the contrary, the neutral element and existence of an opposite are a risk), any commutative semi-group will do; is something more suitable than $\mathbb Z_{2^{512}}$?
- Is $H_i(s_1,H_j(s_2,x))=H_j(s_2,H_i(s_1,x))$ a consequence of $H_1(s_1,H_2(s_2,x)) = H_2(s_2,H_1(s_1,x))$?
- If not, can we make a construction with only the originally asked property?