Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

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

Closed
sirgeorgesawcon opened this issue Aug 29, 2023 · 7 comments
Closed

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

sirgeorgesawcon opened this issue Aug 29, 2023 · 7 comments
Labels
bug Something isn't working

Comments

@sirgeorgesawcon
Copy link

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,

Qiskit/qiskit-ibm-runtime#1043
@sirgeorgesawcon sirgeorgesawcon added the bug Something isn't working label Aug 29, 2023
@ElePT
Copy link
Collaborator

ElePT commented Aug 29, 2023 •
edited
Loading

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
Copy link
Author

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
Copy link
Author

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

@ElePT
Copy link
Collaborator

ElePT commented Aug 30, 2023 •
edited
Loading

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 ElePT mentioned this issue Aug 30, 2023
Closed
@ElePT
Copy link
Collaborator

ElePT commented Aug 30, 2023

See issue here: Qiskit/qiskit#10735

@sirgeorgesawcon
Copy link
Author

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

@ElePT
Copy link
Collaborator

ElePT commented Sep 8, 2023 •
edited
Loading

Issue fixed in Qiskit/qiskit#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.

@ElePT ElePT closed this as completed Sep 8, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
bug Something isn't working
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@ElePT @sirgeorgesawcon