IAE works differently for sampler from Aer and Sampler from qiskit ibm runtime.

sirgeorgesawcon commented 1 year ago

Environment

Qiskit Algorithms version: 0.2.0
Python version: 3.10.12
Operating system: Linux

What is happening?

Describe the bug : I am working on an Iterative Amplitude Estimation routine that works on sampler. I have a quantum circuit, that I run first using the sampler from aer, I import it from there, simply set sampler = _Sampler(run_options={"shots": num_shots})

and run the algorithm, It gives me a result of 0.64165 which is very near to what I shoud be getting (i verified it using classical caluclations , it should be around 0.64)

Then I use from qiskit_ibm_runtime import QiskitRuntimeService, Sampler, Estimator, Options

where my qiskit_ibm_runtime version is 0.11.3

Then I run the exact same quantum circuit qc, using backend = "ibmq_qasm_simulator"

options = Options() options.execution.shots = 2000 options.resilience_level = 0

and then invoke the session and make use of sampler = Sampler(options=options)

The result I'm getting from this is never the same, and is no where close to my answer that I was getting using aer Sampler.

Even though the backend I'm using is IBM_Qasm_simulator , which is error free and ideal simulator, so there should not be any error in the result, but I'm not able to replicate the result that I was getting using the aer sampler. Why is it so? Since both these make use of shot-type backend, the only difference is that one is running on my local device and the other on IBM cloud. But the answers are so different.

Suggested solutions :

Additional Information : I also had a lengthy conversation on Qiskit's slack as well on this said topic Slack Chat: https://qiskit.slack.com/archives/C7SJ0PJ5A/p1693194989748429?thread_ts=1693194989.748429&cid=C7SJ0PJ5A

qiskit-ibm-runtime version: 0.11.3
Python version: 3.10.12
Operating system: Linux

How can we reproduce the issue?

Steps to reproduce : Here's the code for the quantum circuit

n = 3
obj = 8
mu = 4
c = 0.1
sigma = 1
a = 0
b = 7
reps = 2
p_s = 0.5
p_b = 0.2
demand = NormalDistribution(num_qubits=n,mu=mu,sigma=sigma,  bounds=[a,b])
supply = RealAmplitudes(n, reps=reps)
supply = supply.assign_parameters(parameters)
qc = QuantumCircuit(9)
qc.compose(demand,qubits=[3,4,5], inplace=True)
qc.compose(supply, qubits = [0,1,2], inplace=True)
z = 6
z1 = 7
for i in range(n):
    qc.x(i)

for i in range(1,n):
    qc.cnot(z,i)

for i in range(1,n):
    qc.cnot(z,i+n)

qc.ccx(0,n,n+n)

qc.x(0)

for i in range(1,n-1):
  qc.cnot(n+n,i)
  qc.cnot(n+n,i+n)
  qc.ccx(i+n,i,n+n)
  qc.ccx(i+n,i,n)

qc.cnot(n+n,n-1)
qc.cnot(n+n,n+n-1)
qc.ccx(n+n-1,n-1,n+n)

for i in range(1,n-1):
    qc.cnot(n,n-i)
    qc.cnot(n,2*n-i)
    qc.ccx(2*n-i-1,n-i-1,n)

for i in range(1,n-1):
    qc.ccx(0,n,n-i)
    qc.ccx(0,n,2*n-i)

qc.x(0)

qc.ccx(0,n,1)

qc.ccx(0,n,n+1)

for i in range(n):
  qc.x(i)

qc.x(2*n)

for i in range(n):
    qc.x(i+n)

for i in range(1,n):
    qc.cnot(z1,i)

for i in range(1,n):
    qc.cnot(z1,i+n)

qc.ccx(0,n,z1)

qc.x(n)

for i in range(1,n-1):
  qc.cnot(z1,i+n)
  qc.cnot(z1,i)
  qc.ccx(i,i+n,z1)
  qc.ccx(i,i+n,0)

qc.cnot(z1,z1-2)
qc.cnot(z1,n-1)
qc.ccx(z1-2,n-1,z1)

for i in range(1,n-1):
    qc.cnot(0,z1-(i+1))
    qc.cnot(0,n-i)
    qc.ccx(z1-i-2,n-i-1,0)

for i in range(1,n-1):
    qc.ccx(0,n,z1-(i+1))
    qc.ccx(0,n,n-i)

qc.x(n)
qc.ccx(0,n,n+1)
qc.ccx(0,n,1)

for i in range(n):
    qc.x(n+i)
qc.x(z1)

for i in range(n):
    c2ry = RYGate(2*c*p_s * 2**(i)).control(2)
    qc.append(c2ry,[i+n,obj-2,obj])

for i in range(n):
    c2ry = RYGate(2*c*p_s * 2**(i)).control(2)
    qc.append(c2ry,[i,obj-1,obj])

for i in range(n):
    c2ry = RYGate(-2*c*p_s * 2**(i)).control(3)
    qc.append(c2ry,[i,obj-1,obj-2,obj])

for i in range(n):
  qc.cry(-2*c*p_b*(2**i),i,obj)

qc.ry(np.pi/2,obj)
qc.draw('mpl')

Once the quantum circuit qc is made, first we import sampler from :

from qiskit_aer.primitives import Sampler

Then call the IAE function:

epsilon = 0.01   # the error
alpha = 0.5      # the alpha value

problem = EstimationProblem(state_preparation=qc,objective_qubits=[obj])     # the quantum circuit

# construct amplitude estimation
iae = IterativeAmplitudeEstimation(
    epsilon_target=epsilon,
    alpha=alpha,
    sampler=Sampler(run_options={"shots": num_shots})
)

and then:

result = iae.estimate(problem)
print(result)
print((result.estimation-0.5)/c)

The result for the final print command is 0.641

But then, when i make use of :

from qiskit_ibm_runtime import QiskitRuntimeService, Sampler, Estimator, Session, Options
# Save an IBM Quantum account.
QiskitRuntimeService.save_account(channel="ibm_quantum", token="MY API TOKEN",overwrite=True)

and then set the backend:

backend = "ibmq_qasm_simulator"

and then the options

options = Options()
options.execution.shots = num_shots
options.resilience_level = 0

and finally calling the IAE:

# runt he circuit on sampler

epsilon = 0.01   # the error
alpha = 0.5

with Session(service=service, backend=backend):
    sampler = Sampler(options = options)
    iae = IterativeAmplitudeEstimation(epsilon_target=epsilon,alpha=alpha, sampler=sampler)
    problem = EstimationProblem(state_preparation=qc,objective_qubits=[obj])
    result = iae.estimate(problem)
    print(result)
    print((result.estimation-0.5)/c)

always give me a different result from what is Intended.

What should happen?

Expected behavior : I expect to see the same behaviour as I was seeing from the aer sampler in qiskit ibm runtime sampler.

Any suggestions?

I already discussed this issue on Slack : https://qiskit.slack.com/archives/C7SJ0PJ5A/p1693194989748429?thread_ts=1693194989.748429&cid=C7SJ0PJ5A

and also opened an Issue in qiskit-ibm-runtime and they suggested to open an issue here,

https://github.com/Qiskit/qiskit-ibm-runtime/issues/1043

ElePT commented 1 year ago

Hi! I have not been able to reproduce your bug. I filled in your snippets with the necessary imports but the algorithm raises an unexpected error. It would help a lot if you could provide a single complete snippet that I could copy-paste and run to reproduce the behavior you saw (the runtime one is enough, no need for another snippet for the aer sampler).

sirgeorgesawcon commented 1 year ago

sure For the first iteration, when I'm using only the aer sampler:

# doing the basic imports

import numpy as np
import math
import matplotlib.pyplot as plt

# quantum imports

from qiskit import *                    # getting all the necessary qiskit packages

from qiskit import QuantumRegister, QuantumCircuit, ClassicalRegister
from qiskit.visualization import *      # visualization tools
from qiskit.circuit import ParameterVector
from qiskit.circuit.library import IntegerComparator,TwoLocal
from qiskit.circuit.library.arithmetic import LinearAmplitudeFunction, CDKMRippleCarryAdder
from qiskit.circuit.library import RealAmplitudes

simulator = BasicAer.get_backend('qasm_simulator')  # setting the backend
from qiskit_algorithms import IterativeAmplitudeEstimation, EstimationProblem
from qiskit_algorithms.optimizers import optimizer, ADAM,GradientDescent, SPSA, COBYLA, GSLS, SLSQP, AQGD, NELDER_MEAD, POWELL

from qiskit_finance.applications.estimation import EuropeanCallPricing
from qiskit_finance.circuit.library import NormalDistribution
from qiskit.circuit.library import RYGate
from qiskit_aer.primitives import Sampler
from qiskit.quantum_info import Statevector
#from qiskit_aer.primitives import Sampler
from qiskit.utils import QuantumInstance, algorithm_globals

then making the quantum circuit:

n = 3
obj = 8
mu = 4
c = 0.1
sigma = 1
a = 0
b = 7
reps = 2
p_s = 0.5
p_b = 0.2
parameters = [np.pi/2] * n * (reps+1) 
demand = NormalDistribution(num_qubits=n,mu=mu,sigma=sigma,  bounds=[a,b])
supply = RealAmplitudes(n, reps=reps)
supply = supply.assign_parameters(parameters)
qc = QuantumCircuit(9)
qc.compose(demand,qubits=[3,4,5], inplace=True)
qc.compose(supply, qubits = [0,1,2], inplace=True)
z = 6
z1 = 7
for i in range(n):
    qc.x(i)

for i in range(1,n):
    qc.cnot(z,i)

for i in range(1,n):
    qc.cnot(z,i+n)

qc.ccx(0,n,n+n)

qc.x(0)

for i in range(1,n-1):
  qc.cnot(n+n,i)
  qc.cnot(n+n,i+n)
  qc.ccx(i+n,i,n+n)
  qc.ccx(i+n,i,n)

qc.cnot(n+n,n-1)
qc.cnot(n+n,n+n-1)
qc.ccx(n+n-1,n-1,n+n)

for i in range(1,n-1):
    qc.cnot(n,n-i)
    qc.cnot(n,2*n-i)
    qc.ccx(2*n-i-1,n-i-1,n)

for i in range(1,n-1):
    qc.ccx(0,n,n-i)
    qc.ccx(0,n,2*n-i)

qc.x(0)

qc.ccx(0,n,1)

qc.ccx(0,n,n+1)

for i in range(n):
  qc.x(i)

qc.x(2*n)

for i in range(n):
    qc.x(i+n)

for i in range(1,n):
    qc.cnot(z1,i)

for i in range(1,n):
    qc.cnot(z1,i+n)

qc.ccx(0,n,z1)

qc.x(n)

for i in range(1,n-1):
  qc.cnot(z1,i+n)
  qc.cnot(z1,i)
  qc.ccx(i,i+n,z1)
  qc.ccx(i,i+n,0)

qc.cnot(z1,z1-2)
qc.cnot(z1,n-1)
qc.ccx(z1-2,n-1,z1)

for i in range(1,n-1):
    qc.cnot(0,z1-(i+1))
    qc.cnot(0,n-i)
    qc.ccx(z1-i-2,n-i-1,0)

for i in range(1,n-1):
    qc.ccx(0,n,z1-(i+1))
    qc.ccx(0,n,n-i)

qc.x(n)
qc.ccx(0,n,n+1)
qc.ccx(0,n,1)

for i in range(n):
    qc.x(n+i)
qc.x(z1)

for i in range(n):
    c2ry = RYGate(2*c*p_s * 2**(i)).control(2)
    qc.append(c2ry,[i+n,obj-2,obj])

for i in range(n):
    c2ry = RYGate(2*c*p_s * 2**(i)).control(2)
    qc.append(c2ry,[i,obj-1,obj])

for i in range(n):
    c2ry = RYGate(-2*c*p_s * 2**(i)).control(3)
    qc.append(c2ry,[i,obj-1,obj-2,obj])

for i in range(n):
  qc.cry(-2*c*p_b*(2**i),i,obj)

qc.ry(np.pi/2,obj)
qc.draw('mpl')

Then after the ciruit is made :


epsilon = 0.01   # the error
alpha = 0.5      # the alpha value
num_shots = 2000

problem = EstimationProblem(state_preparation=qc,objective_qubits=[obj])     # the quantum circuit

# construct amplitude estimation
iae = IterativeAmplitudeEstimation(
    epsilon_target=epsilon,
    alpha=alpha,
    sampler=Sampler(run_options={"shots": num_shots})
)

and then checking the result:

result = iae.estimate(problem)
print(result)
print((result.estimation-0.5)/c)

It will give out a result of 0.64

Now for the runtime part:

# loading the account

import numpy as np
from qiskit import *
from qiskit_ibm_runtime import QiskitRuntimeService, Sampler, Estimator, Session, Options
# Save an IBM Quantum account.
QiskitRuntimeService.save_account(channel="ibm_quantum", token="API TOKEN HERE",overwrite=True)

# doing the basic imports
import math
import matplotlib.pyplot as plt

from qiskit import QuantumRegister, QuantumCircuit, ClassicalRegister
from qiskit.visualization import *      # visualization tools
from qiskit.circuit import ParameterVector
from qiskit.circuit.library import IntegerComparator,TwoLocal
from qiskit.circuit.library.arithmetic import LinearAmplitudeFunction, CDKMRippleCarryAdder
from qiskit.circuit.library import RealAmplitudes

simulator = BasicAer.get_backend('qasm_simulator')  # setting the backend
from qiskit_algorithms import IterativeAmplitudeEstimation, EstimationProblem
from qiskit_algorithms.optimizers import optimizer, ADAM,GradientDescent, SPSA, COBYLA, GSLS, SLSQP, AQGD, NELDER_MEAD, POWELL

from qiskit_finance.circuit.library import NormalDistribution
from qiskit.circuit.library import RYGate, XGate

from qiskit.quantum_info import Statevector
from qiskit.utils import QuantumInstance, algorithm_globals

from qiskit.transpiler import PassManager, InstructionDurations
from qiskit.transpiler.passes import ALAPSchedule, DynamicalDecoupling
from qiskit.visualization import timeline_drawer
service = QiskitRuntimeService(channel='ibm_quantum')

#service = QiskitRuntimeService()

The same code for the Quantum Circuit i.e

n = 3
obj = 8
mu = 4
c = 0.1
sigma = 1
a = 0
b = 7
reps = 2
p_s = 0.5
p_b = 0.2
demand = NormalDistribution(num_qubits=n,mu=mu,sigma=sigma,  bounds=[a,b])
supply = RealAmplitudes(n, reps=reps)
supply = supply.assign_parameters(parameters)
qc = QuantumCircuit(9)
qc.compose(demand,qubits=[3,4,5], inplace=True)
qc.compose(supply, qubits = [0,1,2], inplace=True)
z = 6
z1 = 7
for i in range(n):
    qc.x(i)

for i in range(1,n):
    qc.cnot(z,i)

for i in range(1,n):
    qc.cnot(z,i+n)

qc.ccx(0,n,n+n)

qc.x(0)

for i in range(1,n-1):
  qc.cnot(n+n,i)
  qc.cnot(n+n,i+n)
  qc.ccx(i+n,i,n+n)
  qc.ccx(i+n,i,n)

qc.cnot(n+n,n-1)
qc.cnot(n+n,n+n-1)
qc.ccx(n+n-1,n-1,n+n)

for i in range(1,n-1):
    qc.cnot(n,n-i)
    qc.cnot(n,2*n-i)
    qc.ccx(2*n-i-1,n-i-1,n)

for i in range(1,n-1):
    qc.ccx(0,n,n-i)
    qc.ccx(0,n,2*n-i)

qc.x(0)

qc.ccx(0,n,1)

qc.ccx(0,n,n+1)

for i in range(n):
  qc.x(i)

qc.x(2*n)

for i in range(n):
    qc.x(i+n)

for i in range(1,n):
    qc.cnot(z1,i)

for i in range(1,n):
    qc.cnot(z1,i+n)

qc.ccx(0,n,z1)

qc.x(n)

for i in range(1,n-1):
  qc.cnot(z1,i+n)
  qc.cnot(z1,i)
  qc.ccx(i,i+n,z1)
  qc.ccx(i,i+n,0)

qc.cnot(z1,z1-2)
qc.cnot(z1,n-1)
qc.ccx(z1-2,n-1,z1)

for i in range(1,n-1):
    qc.cnot(0,z1-(i+1))
    qc.cnot(0,n-i)
    qc.ccx(z1-i-2,n-i-1,0)

for i in range(1,n-1):
    qc.ccx(0,n,z1-(i+1))
    qc.ccx(0,n,n-i)

qc.x(n)
qc.ccx(0,n,n+1)
qc.ccx(0,n,1)

for i in range(n):
    qc.x(n+i)
qc.x(z1)

for i in range(n):
    c2ry = RYGate(2*c*p_s * 2**(i)).control(2)
    qc.append(c2ry,[i+n,obj-2,obj])

for i in range(n):
    c2ry = RYGate(2*c*p_s * 2**(i)).control(2)
    qc.append(c2ry,[i,obj-1,obj])

for i in range(n):
    c2ry = RYGate(-2*c*p_s * 2**(i)).control(3)
    qc.append(c2ry,[i,obj-1,obj-2,obj])

for i in range(n):
  qc.cry(-2*c*p_b*(2**i),i,obj)

qc.ry(np.pi/2,obj)
qc.draw('mpl')

This time the backend is:

backend = "ibmq_qasm_simulator"

and then setting the options:

options = Options()

options.execution.shots = 2000
#options.seed_simulator = 120
options.optimization_level = 0

options.resilience_level = 0

and finally running the circuit:

# runt he circuit on sampler

epsilon = 0.001   # the error
alpha = 0.5

with Session(service=service, backend=backend):
    sampler = Sampler(options = options)
    iae = IterativeAmplitudeEstimation(epsilon_target=epsilon,alpha=alpha, sampler=sampler)
    problem = EstimationProblem(state_preparation=qc,objective_qubits=[obj])
    result = iae.estimate(problem)
    print(result)
    print((result.estimation-0.5)/c)

This will yield a result of around 0.02 , and if you run this again, it'll show something different( like -0057). The answer is so different than the aer sampler, even though the backend is an ideal noise free simulator

sirgeorgesawcon commented 1 year ago

I have also uploaded it here : https://github.com/sirgeorgesawcon/test/blob/main/Copy_of_Error_Issue.ipynb

ElePT commented 1 year ago

Thanks a lot @sirgeorgesawcon. It turns out that the serialization module used to send the information over to the runtime server does not deal well with opaque parameterized instructions such as the multi-controlled RY gates you are using. I will open an issue in qiskit to expose the bug, but in the meantime you can work around the issue if you transpile your circuit before creating your EstimationProblem, this is, running:

from qiskit import transpile

....
# here is where you build your state_preparation circuit
qc = QuantumCircuit(...)
...
for i in range(n):
  qc.cry(-2*c*p_b*(2**i),i,obj)
qc.ry(np.pi/2,obj)
qc.draw('mpl')

# once you are done creating the circuit, transpile
qc = transpile(qc, basis_gates=['sx', 'x', 'rz', 'cx'])

Let me know if this works for you :) I tried it out in my environment and got 0.6394340087783801.

ElePT commented 1 year ago

See issue here: https://github.com/Qiskit/qiskit/issues/10735

sirgeorgesawcon commented 1 year ago

Hi @ElePT , thank you for you help. Yes the code works now.

ElePT commented 1 year ago

Issue fixed in https://github.com/Qiskit/qiskit/pull/10758 :) The change will not be deployed in the runtime environment at least until qiskit's next bugfix release, so it might take some time to see it fixed in the Sampler.

qiskit-community / qiskit-algorithms