Qiskit / qiskit

Qiskit is an open-source SDK for working with quantum computers at the level of extended quantum circuits, operators, and primitives.
https://www.ibm.com/quantum/qiskit
Apache License 2.0
4.85k stars 2.29k forks source link

``qiskit.opflow.PauliTrotterEvolution`` of ``SuzukiTrotter`` has wrong gate-angles for odd orders > 1 #7623

Open fazlekarim-anthem opened 2 years ago

fazlekarim-anthem commented 2 years ago

Environment

What is happening?

When running a PauliTrotterEvolution with odd order greater than 1 on an EvolvedOp operator that has a PauliOp or PauliSumOp primitive, the code ignores the parameters binded to it.

How can we reproduce the issue?

When running the following example:

import qiskit
from qiskit.circuit import Parameter
from qiskit.opflow import X, Y, Zero, PauliTrotterEvolution, Suzuki, CircuitStateFn
from qiskit.circuit.quantumcircuit import QuantumCircuit

statebound_op = CircuitStateFn(QuantumCircuit(1))
theta = Parameter("theta")
ham = theta * (X+Y)
evo_ham = ham.exp_i()
state = Zero ^ Zero
trotter_op = PauliTrotterEvolution(trotter_mode=Suzuki(order=2, reps=1)).convert(evo_ham)
bound_op = trotter_op.bind_parameters({theta: 0.33})
evo_state = bound_op @ statebound_op
evo_state.to_circuit().draw()

where the order is 2 (or an even number) we get the following circuit: image

If we repeat the above code with an order of 3 (or an odd number greater than 2) we will get the following circuit: image

As you notice, the parameters are not being binded when the order is an odd number greater than 2.

What should happen?

The expected outcome is when there should be a parameter binded to the circuit for odd orders

Any suggestions?

(As a disclaimer, I am still familiarizing myself with the subject and codebase).

I came to the realization that regardless of what the trotter_mode is in PauliTrotterEvolution, if the operator is an EvolvedOP and if its primitive is a PauliOp or PauliSumOp, we will use the predefined trotter in the class where its SuzukiTrotter when the order is greater than 1 and LieTrotter when the order is 1 (regardless of what the user defines it to be). This SuzukiTrotter used in the class has a different implementation than qiskit.opflow.evolutions.trotterizations.suzuki.py.

When we look into the implementation of SuzukiTrotter: https://github.com/Qiskit/qiskit-terra/blob/e148f97b9afe869511fe82dc6a0e83b7b33dd217/qiskit/synthesis/evolution/suzuki_trotter.py#L129-L146

we notice when order==1 we return a pauli_list without multiplying it with time. This begs the question, what do we mean by order here? Is an order of n the same as 2k of the implemented paper or an order of n the same as k?

If the answer is the prior, then we need to have an assertion such that the order, n, can not be an odd number as k is always an integer greater than 1.

If its the later, then we must return the following [(op, coeff * time) for op, coeff in pauli_list] when order==1.

I am looking at order == 1 because we recursively will go back to it when the order is an odd number greater than 2: https://github.com/Qiskit/qiskit-terra/blob/e148f97b9afe869511fe82dc6a0e83b7b33dd217/qiskit/synthesis/evolution/suzuki_trotter.py#L139-L146.

One more thing to do is compare these two codes (SuzukiTrotter and qiskit.opflow.evolutions.trotterizations.suzuki.py ). In one implementation time is taken in consideration and in the other time is not.

Cryoris commented 2 years ago

The bound parameter values may not show up in your circuit diagram as being bound, but they are 🙂 Here's your example with 3 repetitions (odd and larger than 2):

   ┌───────────┐┌───────────┐┌───────────┐┌───────────┐┌───────────┐┌───────────┐┌───────────┐┌───────────┐┌───────────┐┌───────────┐
q: ┤ exp(it X) ├┤ exp(it Y) ├┤ exp(it X) ├┤ exp(it Y) ├┤ exp(it X) ├┤ exp(it Y) ├┤ exp(it X) ├┤ exp(it Y) ├┤ exp(it X) ├┤ exp(it Y) ├
   └───────────┘└───────────┘└───────────┘└───────────┘└───────────┘└───────────┘└───────────┘└───────────┘└───────────┘└───────────┘

but if we call .decompose() on the circuit and have a look inside we see the values are actually bound:

   ┌───────┐┌───────┐┌───────┐┌───────┐┌───────┐┌───────┐┌───────┐┌───────┐┌───────┐┌───────┐
q: ┤ Rx(2) ├┤ Ry(2) ├┤ Rx(2) ├┤ Ry(2) ├┤ Rx(2) ├┤ Ry(2) ├┤ Rx(2) ├┤ Ry(2) ├┤ Rx(2) ├┤ Ry(2) ├
   └───────┘└───────┘└───────┘└───────┘└───────┘└───────┘└───────┘└───────┘└───────┘└───────┘

By the way, you might want to consider using the qiskit.circuit.library.PauliEvolutionGate for these evolutions as it's the most up to date implementation 🙂

fazlekarim-anthem commented 2 years ago

^^ In the above example, the values that we set are not bounded. If you see, I am trying to bind the number 0.33, but we are still not seeing its impact. image

I will try using qiskit.circuit.library.PauliEvolutionGate and will get back to you on it.

fazlekarim-anthem commented 2 years ago

It seems like qiskit.circuit.library.PauliEvolutionGate does not allow the operators that have coefficients of type Parameter. It throws an error. I may be using it incorrectly. This is what I have done:

theta = Parameter("theta")
ham = theta * (X+Y)

PauliEvolutionGate(sum_gates)

to get the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/var/folders/dy/qsm_4wdj2dl14m2hpvs42vr00000gn/T/ipykernel_36831/1517703313.py in <module>
      2 ham = theta*(X+Y)
      3 
----> 4 PauliEvolutionGate(sum_gates)

~/PycharmProjects/Anthem-Quantum/qcad/conda/qcad_attempt10/lib/python3.7/site-packages/qiskit/circuit/library/pauli_evolution.py in __init__(self, operator, time, label, synthesis)
     71             name = f"exp(-i {[' + '.join(op.paulis.to_labels()) for op in operator]})"
     72         else:
---> 73             operator = _to_sparse_pauli_op(operator)
     74             name = f"exp(-i {' + '.join(operator.paulis.to_labels())})"
     75 

~/PycharmProjects/Anthem-Quantum/qcad/conda/qcad_attempt10/lib/python3.7/site-packages/qiskit/circuit/library/pauli_evolution.py in _to_sparse_pauli_op(operator)
    104     if isinstance(operator, PauliSumOp):
    105         sparse_pauli = operator.primitive
--> 106         sparse_pauli._coeffs *= operator.coeff
    107         return sparse_pauli
    108     if isinstance(operator, PauliOp):

TypeError: ufunc 'multiply' output (typecode 'O') could not be coerced to provided output parameter (typecode 'D') according to the casting rule ''same_kind''
Cryoris commented 2 years ago

It does support it as time parameter 🙂

theta = Parameter("theta")
ham = X + Y

PauliEvolutionGate(ham, time=theta)
Cryoris commented 2 years ago

Ah now I understand your comment above: the parameters are bound, but to the wrong values!

The Suzuki-Trotter expansion formula should not support odd orders larger than 1, you will find that the PauliEvolutionGate throws a warning if you try that and in the future it will throw an error. That's because, as you also pointed out above, the paper on the Suzuki-Trotter expansion only describes even orders.

fazlekarim-anthem commented 2 years ago

Ahh. I am glad you will throw an error in the future. Just an fyi, I do not see a warning with PauliEvolutionGate with the synthesizer as SuzukiTrotter.

theta = Parameter("theta")
ham = X + Y
op_exp = PauliEvolutionGate(ham, time=theta, synthesis=SuzukiTrotter(order= 3, reps= 1))
circ = QuantumCircuit(1)
circ.append(PEG,range(1))
print(circ.decompose())

leads to the following output again without any warning: image

Further, it seems like the time input can take in one parameter (theta). If I wanted alphaX+betaY (or the coefficient of X to be alpha and coefficient of Y to be beta), I would not be able to use PauliEvolutionGate. What would I do then?

Cryoris commented 2 years ago

Ah, the warning will only come with the dev version and the 0.20 release of Terra 🙂

Further, it seems like the time input can take in one parameter (theta). If I wanted alphaX+betaY (or the coefficient of X to be alpha and coefficient of Y to be beta), I would not be able to use PauliEvolutionGate. What would I do then?

That's currently not supported you're right. That would require SparsePauliOps with coefficients of type ParameterExpressions -- luckily this is coming with #7215. After that's merged we should be able to cover your use case.