Speed up tr(H^k) by another k! factor #80

Context

Computing high order moments $\mu_k=\mathrm{tr}(H^k)$ is important in statistical physics, for example to compute high temperature expansions: all the information about the thermodynamic properties of a Hamiltonian $H$ is contained in the moments $\mu_k$.

The Pauli strings representation is quite efficient are computing these because it exploits the fact that $H$ is sparse in the Pauli basis, and we do not need to construct the full $H^k$ in order to access $\mathrm{tr}(H^k)$.
The current implementation is based on decomposing $\mathrm{tr}(H^k)=\mathrm{tr}(H^{k/2}H^{k/2} )$ and computing the norm of $H^{k/2}$.
It seems however that there is a more efficient way to do this:
Expanding $\mathrm{tr}(H^k) = \sum_{l_1, \dots, l_k} c_{l_1} \cdots c_{l_k} , \mathrm{tr}(P_{l_1} \cdots P_{l_k})$, the trace inside is either $0$ or in ${\pm 1, \pm i}$. Whether it vanishes depends only on the multiset (the same Pauli can appear multiple times) of Paulis, not their order. If one ordering gives zero, all $k!$ orderings give zero.
(cf. #51 (comment))

Goal

Exploit this by (1) enumerating only the multisets ${P_{l_1}, \dots, P_{l_k}}$ whose product is proportional to the identity, (2) computing the phase analytically for each multiset, (3) summing the contributions. This should save a factor of $k!$ over the naive sum.

• Write an algorithm to enumerate the identity-yielding multisets
• For a fixed multiset, design an efficient way (ideally a closed form formula) to compute $\sum_\sigma \mathrm{phase}(P_{l_{\sigma(1)}} \cdots P_{l_{\sigma(k)}})$ over permutations $\sigma$
• Compare the performance to the existing trace_product

A valid solution should outperfom the current example for the example below and higher order moments.

Application example

Computing the moments of the transverse field Ising model

using PauliStrings

function ising(h, N)
    H = Operator(N)
    for i in 1:N
        H += h, "X", i
    end
    for i in 1:N
        H += "Z", i, "Z", mod1(i + 1, N)
    end
    return H
end

N = 20
H = ising(0.5, N)
for k in 1:14
    muk = trace_product(H,k; scale=1)
    println("mu_", k, " = ", muk)
end