cirq.PauliStringPhasor returns wrong unitary

[pgoiporia]

Description of the issue
I am trying to check the unitary of acirq.PauliStringPhasor. I don't get what I expect.

How to reproduce the issue

import cirq

qubits = cirq.LineQubit.range(2)
phasor = cirq.PauliStringPhasor(cirq.Z(qubits[0]), qubits, exponent_neg=1)
circuit = cirq.Circuit(phasor)

print(circuit)
print(cirq.unitary(circuit))

The unitary this generates represents a Z⊗Z, not Z⊗I

Cirq version
1.4.0.dev20240403011731

[burlemarxiste]

I stumbled upon this when investigating a related request.
#6598

Here is a proposed fix, with a unit test that checks the unitaries for all "dense" Pauli strings of len <= 4.
burlemarxiste@09d37f4

As for why this extra Z...

The algorithm used in PauliStringPhasorGate._decompose_ is detailed in section 4.7.3 of N&C (with the only difference that it is done in place instead of on an ancilia):

• We perform a change of basis in which all X and Y in the Pauli string are substituted by Z.
• A CNOT computes the parity of each basis state.
• If it's even the state has to be phase-shifted by exponent_pos, otherwise by exponent_neg
• We undo the CNOT and the change of basis.

The bug was caused by wrongly involving the extra qubit in the CNOT. The proposed fix partitions qubits into a set affected by this process (those from the Pauli string) and another left aside. A sequence of I gates are applied to the other qubits, for the sole purpose of displaying a unitary of the correct size (otherwise the tensor products by $I_2$ are implicit).

Finally, in case the argument is (or simplifies to) I, a global phase shift is still applied. This might seem superfluous, but this ensures a correct behavior when this gate is controlled (taking advantage of the fact that GlobalPhase gates, when controlled;, are turned into a controlled phase).

[NoureldinYosri]

oh I forgot to reply to this... it seems what you are looking for is cirq.DensePauliString not cirq. PauliStringPhasor. cirq.PauliStringPhasor adds extra phases see https://quantumai.google/reference/python/cirq/PauliStringPhasor

[burlemarxiste]

I'm afraid this is still an issue!

cirq.PauliStringPhasor(a_pauli_string, qubits, exponent_neg=1)

Should phase the positive eigenvalues (default argument exponent_pos=1) and negative eigenvalues (exponent_neg=1 here) by the same amount, ie, it only applies a global phase shift of $e^i\pi$, so the result should have the same unitary as the a_pauli_string argument (up to a global phase shift).

EDIT: Should leave the positive eigenvalues untouched (default argument exponent_pos=0) and set the negative eigenvalues (exponent_neg=1 here) to -1, so the result should have the same unitary as the a_pauli_string argument (up to a global phase shift).

This is probably why @pgoiporia chose these arguments to make the issue obvious to see.

Up to a global phase, the expected unitary for cirq.PauliStringPhasor(cirq.Z(qubits[0]), qubits, exponent_neg=1) should thus be $Z \otimes I$.

The code returns $Z \otimes Z$ here because of the issue I have mentioned causing extra qubits to be incorrectly involved in the parity computation (the CNOTs), and thus be treated as acted upon an implicit Z.

@pgoiporia, do you have other test cases you want to try?

[NoureldinYosri]

I agree that the unitary is wrong but the expected unitary for exponent_neg=1 is not $Z \otimes I$ but $I \otimes I$. the default argument for exponent_pos is zero (not 1)

Cirq/cirq-core/cirq/ops/pauli_string_phasor.py

Line 70 in ab96766

exponent_pos: Union[int, float, sympy.Expr] = 0,

so nothing happens there. for the negative eigen vectors we need to multiply them by -1 which essentially turns the $Z$ into an $I$.

The unitary of this gate is computed from its decomposition which seems to be incorrect.

Cirq/cirq-core/cirq/ops/pauli_string_phasor.py

Lines 354 to 371 in ab96766

	def _decompose_(self, qubits: Sequence['cirq.Qid']) -> 'cirq.OP_TREE':
	if len(self.dense_pauli_string) <= 0:
	return
	any_qubit = qubits[0]
	to_z_ops = op_tree.freeze_op_tree(self._to_z_basis_ops(qubits))
	xor_decomp = tuple(xor_nonlocal_decompose(qubits, any_qubit))
	yield to_z_ops
	yield xor_decomp
	if self.exponent_neg:
	yield pauli_gates.Z(any_qubit) ** self.exponent_neg
	if self.exponent_pos:
	yield pauli_gates.X(any_qubit)
	yield pauli_gates.Z(any_qubit) ** self.exponent_pos
	yield pauli_gates.X(any_qubit)
	yield protocols.inverse(xor_decomp)
	yield protocols.inverse(to_z_ops)

[burlemarxiste]

I double checked.

It is designed to map $Z$ to $diag(e^{i \pi p}, e^{i \pi m})$, not to $diag(+e^{i \pi p}, -e^{i \pi m})$ (as the wording "phasing of the eigenvalues" would suggest).

With $p = 0$ and $m = 1$, the expected result is $Z$ not $I$.

The decomposition is correct and a textbook result, the only problem in the issue raised here is involving the extra qubits in it.

Note that with $p = -m = \theta$ it is effectively turning $Z$ into $diag(e^{i \pi \theta}, e^{-i \pi \theta}) = e^{i \pi \theta Z}$, hence its use for exponentials/rotations.

For a matrix with identical eigenvalues like $I$, one could argue that it should return $diag(e^{i \pi \lambda p}, e^{i \pi \lambda p})$. It then gives consistent results for exponentiation.

[NoureldinYosri]

I added the triage/discuss label to talk about this in the next cirq-sync ... I though the definition of the gate applies the pauli string and then the extra phase .. but per @burlemarxiste it replaces the eigenvalues of the original paulistring but keeps the eigenvectors.. eitherway something incorrect is happening.

[burlemarxiste]

I though the definition of the gate applies the pauli string and then the extra phase

No, the gate evaluates "exponential-like" expressions (exponentials when the positive and negative exponents are opposite).

The code implements the circuit p. 210 of Nielsen&Chuang, with two differences:

• Not using an ancilia: we can compute the parity "in place" on any of the qubits.
• Different phases for the positive and negative eigenspace.

From this circuit it is clear that involving a qubit in the CNOT is as if it had a Z in the pauli string.

burlemarxiste@09d37f4 does everything necessary to make it behave correctly with the identities. The unit test checks the unitaries against values computed by scipy.linalg.expm.

--

I hope you understand now why I focused my attention on this class: it is designed to be a gate to exponentiate Pauli strings, not to apply a Pauli String and then an extra phase, and from the Cirq documentation, it appears to be the "canonical" way of exponentiating a string, not PauliSumExponential:

https://quantumai.google/cirq/build/pauli_observables#exponentials_of_pauli_operators_as_cirqpaulistringphasors

Since PauliStringPhasor has since been updated to allow extra qubits outside the pauli strings with an argument, we might as well make it work as expected with those.

[NoureldinYosri]

@burlemarxiste feel free to open a PR to fix the PauliStringPhasor issue though this won't count towards the fix for the other issue. The other issue needs to be fixed using DensePauliString

this code block won't be valid after we change DensePauliString

        acted_upon_qubits = self.dense_pauli_string.on(*qubits).qubits
        acted_upon_qubit_set = set(acted_upon_qubits)
        identity_qubit_set = set(qubits) - acted_upon_qubit_set

instead you may use an ancilla qubit (same as in the textbook). for that replace _decompose_ with _decompose_with_context_ as

    def _decompose_with_context_(
        self, qubits: Sequence[cirq.Qid], context: Optional[cirq.DecompositionContext] = None
    ) -> cirq.OP_TREE:
        if context is None:
            context = cirq.DecompositionContext(cirq.ops.SimpleQubitManager())
        (ancilla,) = context.qubit_manager.qalloc(1)
        ... # decomposition
       context.qubit_manager.qfree([ancilla])
       
        

[burlemarxiste]

I don't see why the ancilia is needed.

this code block won't be valid after we change DensePauliString

acted_upon_qubits = self.dense_pauli_string.on(*qubits).qubits

Agree, but in the current codebase it is how we get the set of "non-I" qubits from the DensePauliString.

In case of a serious overhaul of the PauliString codebase, I agree that there should be methods to split qubits into a "trivial" and "non-trivial" set. This would be a perfect use case for them!

I'm all for a PauliString overhaul. With a proper design review. Maybe the cirq team could consider making it a Summer of Code project! 😉

[burlemarxiste]

One of these:

A.

The expected behavior of PauliStringPhasor is to preserve the eigenvectors and replace the eigenvalues with a value based solely on the parity of the state's bitstring, thus to map $I$ to $diag(e^{i \pi p}, e^{i \pi m})$, exactly the same way it affects $Z$.

In that case:

• Current implementation is correct, $I$ should be handled like $Z$.
• A single PauliStringPhasor cannot implement exponentiation of a PauliString mixing $XYZ$ and $I$. The exponentiation has to be decomposed into a PauliStringPhasor with $p = \theta, m = -\theta$ applied to the $XYZ$ in the string, and a PauliStringPhasor with $p = \theta, m = \theta$ applied to the remaining $I$ qubits.

B.

The expected behavior of PauliStringPhasor is to preserve the eigenvectors and replace each eigenvalue in the tensor product by $e^{i \pi p}$ if it's equal to $1$, $e^{i \pi m}$ if it's $-1$. It will thus map $I$ to $diag(e^{i \pi p}, e^{i \pi p})$ because $I$ has identical eigenvalues, and $Z$ to $diag(e^{i \pi p}, e^{i \pi m})$.

In that case:

• @pgoiporia is right to report an issue, the current implementation incorrectly affects the $I$ space, and it needs to include a global phase shift whenever the input argument is only $I$s (the change I suggested).
• A single PauliStringPhasor, once fixed, will correctly implement exponentiation of a PauliString mixing $XYZ$ and $I$.

C.

The behavior of PauliStringPhasor is considered undefined on $I$.

In that case:

• Current implementation should not have the extra qubits arguments and should check that the string it receives only contains XYZ.
• Exponentiating an arbitrary PauliString has to be decomposed into a PauliStringPhasor on the XYZ part, identities on the rest, and if the XYZ part is empty, a global phase shift.

---

This is particularly relevant if PauliString will support $I$ in the longer term, because expressions like np.exp(1j * cirq.X(q0) * cirq.Y(q1)) rely on a single call to PauliStringPhasor.

This is where the exponentiation happens:

Cirq/cirq-core/cirq/ops/pauli_string.py

Line 836 in 543d9cd

def __rpow__(self, base):

>>> cirq.unitary(np.exp(1j * cirq.X(q0) * cirq.X(q1) * np.pi / 4))
array([[0.70710678+0.j        , 0.        +0.j        ,
        0.        +0.j        , 0.        +0.70710678j],
       [0.        +0.j        , 0.70710678+0.j        ,
        0.        +0.70710678j, 0.        +0.j        ],
       [0.        +0.j        , 0.        +0.70710678j,
        0.70710678+0.j        , 0.        +0.j        ],
       [0.        +0.70710678j, 0.        +0.j        ,
        0.        +0.j        , 0.70710678+0.j        ]])
>>> cirq.unitary(np.exp(1j * cirq.X(q0) * cirq.X(q0) * np.pi / 4))
array([[1.+0.j]])
>>> 

Getting the correct result for the second expression will involve more than changing the storage format of PauliString to handle I: it'll need one of the three options I have mentioned.