LieTrotter does not seem to preserve the gate order even with preserve_order=True

[t-imamichi]

Environment

• Qiskit version: 1.3.2, main
• Python version: 3.12.8
• Operating system: macOS 15.2

What is happening?

When I generate PauliEvolutionGate of XI+ZZ , it corresponds to U(t)=exp(-it (XI+ZZ)) and is approximated as exp(-it XI) exp(-it ZZ) because LieTrotter with reps=1 and preserve_order=True is applied as default.
I think a circuit with the gate represents U(t)|0> = exp(-it XI) exp(-it ZZ) |0> and ZZ part will be applied earlier than XI part.
But Qiskit generates RZZ gate first as follows. The reverse order matches the numerically calculated matrix.

Reference: LieTrotter documentation

How can we reproduce the issue?

from qiskit import QuantumCircuit, __version__
from qiskit.circuit import Parameter
from qiskit.circuit.library import PauliEvolutionGate
from qiskit.quantum_info import Operator, SparsePauliOp
from scipy.linalg import expm

print(__version__)
qc = QuantumCircuit(2)
t = Parameter("t")
op = SparsePauliOp(["XI", "ZZ"])
print(op)
# exp(-it (XI + ZZ))
evo = PauliEvolutionGate(op, t)
qc.append(evo, qargs=qc.qubits)
print("original order")
print(qc)
# exp(-it (XI + ZZ)) -> expect exp(-it XI) exp(-it ZZ)
print(qc.decompose())

# reverse order
qc2 = QuantumCircuit(2)
evo = PauliEvolutionGate(op[::-1], t)
qc2.append(evo, qargs=qc2.qubits)
print("reverse order")
print(qc2.decompose())

# t = 1
reference = expm(-1j * op[0].to_matrix()) @ expm(-1j * op[1].to_matrix())
evo_op = Operator(qc.assign_parameters([1]))
evo_op2 = Operator(qc2.assign_parameters([1]))
print(f"original order: equiv to ref? {evo_op.equiv(reference)}")
print(f"reverse order: equiv to ref? {evo_op2.equiv(reference)}")

output

2.0.0.dev0+315d6ad
SparsePauliOp(['XI', 'ZZ'],
              coeffs=[1.+0.j, 1.+0.j])
original order
     ┌────────────────────────┐
q_0: ┤0                       ├
     │  exp(-it (XI + ZZ))(t) │
q_1: ┤1                       ├
     └────────────────────────┘

q_0: ──────────────■──────────
     ┌───────────┐ │ZZ(2.0*t)
q_1: ┤ Rx(2.0*t) ├─■──────────
     └───────────┘
reverse order

q_0: ─■───────────────────────
      │ZZ(2.0*t) ┌───────────┐
q_1: ─■──────────┤ Rx(2.0*t) ├
                 └───────────┘
original order: equiv to ref? False
reverse order: equiv to ref? True

What should happen?

I expect to apply RZZ earlier than RX if preserve_order=True

Any suggestions?

It might be resolved by reversing the input operator as I did in the script.

[jakelishman]

I don't have an opinion on which way is more expected, and which way Qiskit should output. I just wanted to add: I can easily see arguments for either interpretation, so we should make very clear in the documentation which it is that we mean.

For example:

• I write op = SparsePauliOp(["XI", "ZZ"]), so I would write $\exp\bigl(it (XI + ZZ)\bigr) \approx \exp(itXI)\exp(itZZ)$, and so I expect that the $ZZ$ term is applied to the state first, and then $XI$.
• I wrote op = SparsePauliOp(["XI", "ZZ"]), and I expect the terms to be applied in the order of for term in op: apply(term), so I expect $XI$ be applied first, and then $ZZ$.

I think both are reasonable viewpoints a person could hold, so we should be really super clear what we mean in the documentation.

[Cryoris]

Agreed, though personally I'd probably think of the first one since I would associate the SparsePauliOp with a sum. There's even proof that I did think this, in the following docstring 😬

qiskit/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py

Line 806 in d03a61e

H = Z_1 X_4 + 2 Y_0 Y_3

From a practical perspective, however, changing the order in the PauliEvolutionGate is a bit annoying as it will change everybody's output circuits.

[jakelishman]

Oh yeah, I think both are totally valid. If pressed, I think I'd say that I expect my second option marginally more than my first, which (iiuc) agrees with Qiskit's current behaviour.

I've no strong horse in this race (so I've nothing more to say, I think), other than my weak expectation and general desire for stability. I just wanted to comment that whichever way you choose to resolve this, I think we need very very clear documentation about it.

[t-imamichi]

I agree that changing the current behavior is annoying. It would be good to (1) write a note about the operator order in the documentation and (2) change the draw of PauliEvolutionGate to be aligned with the actual operator order.

As for (2), Qiskit draws the circuit as follows.

op = SparsePauliOp(["XI", "ZZ"])
evo = PauliEvolutionGate(op, t)
qc.append(evo, qargs=qc.qubits)
print(qc)

# output
     ┌────────────────────────┐
q_0: ┤0                       ├
     │  exp(-it (XI + ZZ))(t) │
q_1: ┤1                       ├
     └────────────────────────┘

That's why I expected $$\exp(-it (XI + ZZ)) \approx \exp(-it XI) \exp(-it ZZ)$$. It would be clear to display exp(-it (ZZ + XI))(t) instead.

[Cryoris]

It's the eternal dilemma of stability vs. changing to what we want it to be 😅 changing is annoying but writing yet another note of "Qiskit does something reversed to what you think it is" is also not very attractive...