# Is there an efficient algorithm for decomposing an arbitrary Hamiltonian into Pauli strings?

Basically the title. If I have a $$2^N\times 2^N$$ Hamiltonian $$H$$ of random numbers (we can take the Hamiltonian as normalized if we want) and $$N$$ is an integer, is there an efficient way of writing $$H = \sum_{i}{\beta_iP_i}$$ where $$\beta_i \in \mathbb{C}$$ and $$P_i$$ is a Pauli string and the sum ranges over all such possible tensor product combinations of the Pauli group, $$\{I,X,Y,Z\}$$. Apologies if my notation is non-standard; let me know if any clarification is needed.

Also, I say efficient because I am aware of such solutions as this and this, but these seem to become exponentially hard in $$N$$.

• Certainly, this seems impossible unless there is some additional structure imposed on $H$. Asking for $H$ in the Pauli string basis, given an $H$ in the standard basis requires matrix multiplication, which is why you will always be polynomial in the dimension of the matrix (i.e. exponential in $N$). Commented Oct 3, 2022 at 19:28
• @Condo While it is true that it is exponentially hard -- just specifying H already requires an exponential number of parameters, it does not require matrix multiplication (at least not in the sense of multiplying two $2^N\times 2^N$ matrices). Commented Oct 3, 2022 at 22:03
• The most efficient way (if $H$ is just a full matrix) is to vectorize $H$ (which makes it a $4^N$-component vector) and then apply a basis transformation to the Pauli basis on each of the $4$-dimensional blocks; this can be achieved as a sequence of multiplication of a $4\times 4^{N-1}$ matrix with a $4\times 4$ matrix. Commented Oct 3, 2022 at 22:04
• @NorbertSchuch Ah yes, I completely mistyped what I meant. What I was trying to say was that the change of basis transformation requires matrix multiplication. Commented Oct 5, 2022 at 20:51

Your output size is $$4^n$$ numbers, so you're certainly not going to do better than $$\Omega(4^n)$$ complexity. At least, not without some sort of sparsity guarantee and some way to compute those big terms analogous to the sparse fast fourier transform.
Anyways, it turns out that you can get all $$4^n$$ numbers by interpreting your matrix as a vector and simulating applying the same quantum circuit you would apply to measure in the Bell basis. This takes $$O(4^n)$$ space and $$O(n 4^n)$$ time:
In more classical computing terms: take your matrix of coefficients, permute each column so that the entry at offset $$j$$ in column $$k$$ ends up at offset $$j \oplus k$$ in column $$k$$, then apply the classical Hadamard transform to each column (or maybe each row? I can never keep that straight). The matrix now contains the Pauli term coefficients.
For example see the one proposed recently in https://arxiv.org/abs/2401.16378. In my tests, I could decompose random $$N=11$$-qubit Hamiltonian in about $$\sim 150$$ minutes.