Odd behavior with Qiskit Pauli operators

Question

I was creating a simple implementation of the Hadamard test when I came across the following, what seems like, strange behavior.

Consider the following snippet, which is part of the computating the expected value of the Pauli-$X$ operator with respect to the state $|\psi\rangle = \frac{1}{2}|0\rangle + \frac{\sqrt{3}}{2}|1\rangle$. After measuring the ancillary, we compute $Prob(|0\rangle) - Prob(|1\rangle) = Re(\langle \psi |X|\psi\rangle)$ (see https://en.wikipedia.org/wiki/Hadamard_test_(quantum_computation))

Circuit 1:

reg = QuantumRegister(2)
classical = ClassicalRegister(1)
qc = QuantumCircuit(reg, classical)

def op():
    r = QuantumRegister(1)
    gate = QuantumCircuit(r)
    gate.append(Operator(Pauli(label='X')), [*r])
    return gate.to_gate()

qc.initialize([1/2, np.sqrt(3)/2], [reg[1]])
qc.barrier()
qc.h(reg[0])
qc.append(op().control(),[*reg])
qc.h(reg[0])
qc.barrier()
qc.measure(reg[0], classical)
job = execute(qc, Aer.get_backend('qasm_simulator'), shots=1024)
print(job.result().get_counts())

The result from QASM is $\{'0': 499, '1': 525\}$, which is not correct. I verified with Qiskit's snapshot expectation that expected value is $\frac{\sqrt{3}}{2}$. This was also verified through density matrix computation with numpy.

However, if I replace this with the $CX$ as opposed to the Pauli-$X$ from the quantum info package of qiskit.

Circuit 2

reg = QuantumRegister(2)
classical = ClassicalRegister(1)
qc = QuantumCircuit(reg, classical)

qc.initialize([1/2, np.sqrt(3)/2], [reg[1]])
qc.barrier()
qc.h(reg[0])
qc.cx(reg[0], reg[1])
qc.h(reg[0])
qc.barrier()
qc.measure(reg[0], classical)
job = execute(qc, Aer.get_backend('qasm_simulator'), shots=1024)
print(job.result().get_counts())

I got the following: $\{'0': 962, '1': 62\}$, which taking $(962-62)/1024$ produces an approximation to the correct result.

Next I printed out the transpilation of Circui1

t = transpile(qc, Aer.get_backend('qasm_simulator'), basis_gates=['u1','u2','u3','cx'])
t.draw(output='mpl')

I verified that this circuit does in fact produce the same result as Circuit 1. However, I then did a density matrix computation using numpy, which also produced the same result. However, this revealed that there is an incorrect phase difference, not a global phase difference. If we view the bipartite system as $reg[0] \otimes reg[1]$ where $reg[0]$ is the ancillary control bit, we get a unitary matrix for the gates in between the second set of barriers above (not including the Hadamards) like the following: $\begin{bmatrix} -1 & 0 &0 & 0 \\ 0 & -1 &0 & 0 \\ 0 & 0 &0 & -i \\ 0 & 0 &-i & 0 \\ \end{bmatrix}$

I verified that my computation in numpy was correct through multiple test; as a sanity check I compared my the matrices in the numpy simulation against those printed out by qiskit in the transpilation. So this circuit is not the same as $CX$.

Any ideas as to why I am seeing this?

I thought this might be a Qiskit issue, but I posted it here in case I am missing something, which most likely the case.

Answer 1

It is an issue. Adding control() to the gate introduces a phase difference. You can verify that:

from qiskit.quantum_info import Operator, Pauli

gate = QuantumCircuit(1)
gate.append(Operator(Pauli(label='X')), [0])
gate = gate.control() 

print(Operator(gate).data)

----
Output:

[[ 1   0   0   0]
 [ 0   0   0  1j]
 [ 0   0   1   0]
 [ 0  1j   0   0]]

So circuit 1 actually computes $ \text{Re}\big(\text{i} \cdot\langle \psi| X |\psi \rangle \big) = 0 $