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

[qiskit-community]

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

[sirgeorgesawcon]

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

[ElePT]

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]

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]

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

[ElePT]

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]

See issue here: Qiskit/qiskit#10735

[ElePT]

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.