# How can I Combine Qiskit's VQE Tutorial and TSP Simulator to solve the TSP on a real backend efficiently?

I have a university project where I am trying to solve the Traveling Salesman Problem by using a real quantum backend, rather than just the VQE simulator, as in this tutorial.

I found this code on Qiskit's documentation, about the VQE, and using it with a real backend. I decided to try to amalgamate the two, changing the Hamiltonian in the real backend code, to be the Hamiltonian from the TSP simulation code. I have been waiting for about 4 hours and gotten around 150 iterations and 150 different ansatz energies. I would like to change my x0 value for the start of the iterations to be more optimal, but I don't really understand what the input values to the ansatz even mean, and I also don't understand how to compare the energy output (which is negative for some reason?) to a brute force Classical technique, where an output would be the shortest distance of the TSP.

Please can you help me understand how to make this run more efficiently so I can actually implement my code, and also to understand better how the ansatz parameters can be tweaked? I would like to make use of the library tsp, but it has been deprecated. I kept the code in my repo as similar to the tutorials as possible to avoid errors, and left the Max-cut section. The relevant code is in the last, or bottom block.

from qiskit.circuit.library import EfficientSU2

# ## Step 1: Map classical inputs to a quantum problem
#
# Here we define the problem instance for our VQE algorithm. Although the problem in question can come from a variety of domains, the form for execution through Qiskit Runtime is the same. Qiskit provides a convenience class for expressing Hamiltonians in Pauli form, and a collection of widely used ansatz circuits in the [qiskit.circuit.library](https://docs.quantum-computing.ibm.com/api/qiskit/circuit_library).
#
# Here, our example Hamiltonian is derived from a quantum chemistry problem
#

# In[4]:

hamiltonian = qubitOp

# Our choice of ansatz is the EfficientSU2 that, by default, linearly entangles qubits, making it ideal for quantum hardware with limited connectivity.
#

# In[5]:

ansatz = EfficientSU2(hamiltonian.num_qubits)
ansatz.decompose().draw("mpl", style="iqp")

# From the previous figure we see that our ansatz circuit is defined by a vector of parameters, $$\theta_{i}$$, with the total number given by:
#

# In[6]:

num_params = ansatz.num_parameters
num_params

# ## Step 2: Optimize problem for quantum execution.
#
# We can schedule a series of [qiskit.transpiler](https://docs.quantum-computing.ibm.com/api/qiskit/transpiler) passes to optimize our circuit for a selected backend. This includes a few components:
#
# *   [optimization_level](https://docs.quantum-computing.ibm.com/api/qiskit/transpiler_preset#preset-pass-manager-generation): The lowest optimization level just does the bare minimum needed to get the circuit running on the device; it maps the circuit qubits to the device qubits and adds swap gates to allow all 2-qubit operations. The highest optimization level is much smarter and uses lots of tricks to reduce the overall gate count. Since multi-qubit gates have high error rates and qubits decohere over time, the shorter circuits should give better results.
# *   [Dynamical Decoupling](https://docs.quantum-computing.ibm.com/api/qiskit-ibm-runtime/qiskit_ibm_provider.transpiler.passes.scheduling.PadDynamicalDecoupling): We can apply a sequence of gates to idling qubits. This cancels out some unwanted interactions with the environment.
#

# In[7]:

from qiskit.transpiler import PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.transpiler.passes import (
ALAPScheduleAnalysis,
ConstrainedReschedule,
)
from qiskit.circuit.library import XGate

target = backend.target
pm = generate_preset_pass_manager(target=target, optimization_level=3)
pm.scheduling = PassManager(
[
ALAPScheduleAnalysis(target=target),
ConstrainedReschedule(target.acquire_alignment, target.pulse_alignment),
target=target, dd_sequence=[XGate(), XGate()], pulse_alignment=target.pulse_alignment
),
]
)

ansatz_ibm = pm.run(ansatz)

# In[8]:

ansatz_ibm.draw(output="mpl", idle_wires=False, style="iqp")

# We can also use apply_layout to transform our virtual observable to physical
#

# In[9]:

hamiltonian_ibm = hamiltonian.apply_layout(ansatz_ibm.layout)
hamiltonian_ibm

# ## Step 3: Execute using Qiskit Primitives.
#
# Like many classical optimization problems, the solution to a VQE problem can be formulated as minimization of a scalar cost function.  By definition, VQE looks to find the ground state solution to a Hamiltonian by optimizing the ansatz circuit parameters to minimize the expectation value (energy) of the Hamiltonian.  With the Qiskit Runtime [Estimator](https://docs.quantum-computing.ibm.com/api/qiskit-ibm-runtime/qiskit_ibm_runtime.Estimator) directly taking a Hamiltonian and parameterized ansatz, and returning the necessary energy, the cost function for a VQE instance is quite simple:
#

# In[10]:

def cost_func(params, ansatz, hamiltonian, estimator):
"""Return estimate of energy from estimator

Parameters:
params (ndarray): Array of ansatz parameters
ansatz (QuantumCircuit): Parameterized ansatz circuit
hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
estimator (Estimator): Estimator primitive instance

Returns:
float: Energy estimate
"""
energy = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
return energy

# Note that, in addition to the array of optimization parameters that must be the first argument, we use additional arguments to pass the terms needed in the cost function.
#

# ### Creating a callback function
#
# Callback functions are a standard way for users to obtain additional information about the status of an iterative algorithm.  The standard SciPy callback routine allows for returning only the interim vector at each iteration.  However, it is possible to do much more than this.  Here, we show how to use a mutable object, such as a dictionary, to store the current vector at each iteration, for example in case we need to restart the routine due to failure, and also return the current iteration number and average time per iteration.
#

# In[11]:

def build_callback(ansatz, hamiltonian, estimator, callback_dict):
"""Return callback function that uses Estimator instance,
and stores intermediate values into a dictionary.

Parameters:
ansatz (QuantumCircuit): Parameterized ansatz circuit
hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
estimator (Estimator): Estimator primitive instance
callback_dict (dict): Mutable dict for storing values

Returns:
Callable: Callback function object
"""

def callback(current_vector):
"""Callback function storing previous solution vector,
computing the intermediate cost value, and displaying number
of completed iterations and average time per iteration.

Values are stored in pre-defined 'callback_dict' dictionary.

Parameters:
current_vector (ndarray): Current vector of parameters
returned by optimizer
"""
# Keep track of the number of iterations
callback_dict["iters"] += 1
# Set the prev_vector to the latest one
callback_dict["prev_vector"] = current_vector
# Compute the value of the cost function at the current vector
current_cost = (
estimator.run(ansatz, hamiltonian, parameter_values=current_vector).result().values[0]
)
callback_dict["cost_history"].append(current_cost)
# Print to screen on single line
print(
"Iters. done: {} [Current cost: {}]".format(callback_dict["iters"], current_cost),
end="\r",
flush=True,
)

return callback

# In[12]:

callback_dict = {
"prev_vector": None,
"iters": 0,
"cost_history": [],
}

# We can now use a classical optimizer of our choice to minimize the cost function. Here, we use the [COBYLA routine from SciPy through the minimize function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html). Note that when running on real quantum hardware, the choice of optimizer is important, as not all optimizers handle noisy cost function landscapes equally well.
#
# To begin the routine, we specify a random initial set of parameters:
#

# In[13]:

x0 = 2 * np.pi * np.random.random(num_params)

# Because we are sending a large number of jobs that we would like to execute together, we use a [Session](https://docs.quantum-computing.ibm.com/api/qiskit-ibm-runtime/qiskit_ibm_runtime.Session) to execute all the generated circuits in one block.  Here args is the standard SciPy way to supply the additional parameters needed by the cost function.
#

# In[14]:

# To run on local simulator:
#   1. Use the Estimator from qiskit.primitives instead.
#   2. Remove the Session context manager below.

options = Options()
options.transpilation.skip_transpilation = True
options.execution.shots = 10000

with Session(backend=backend):
estimator = Estimator(options=options)
callback = build_callback(ansatz_ibm, hamiltonian_ibm, estimator, callback_dict)
res = minimize(
cost_func,
x0,
args=(ansatz_ibm, hamiltonian_ibm, estimator),
method="cobyla",
callback=callback,
)

# At the terminus of this routine we have a result in the standard SciPy OptimizeResult format.  From this we see that it took nfev number of cost function evaluations to obtain the solution vector of parameter angles (x) that, when plugged into the ansatz circuit, yield the approximate ground state solution we were looking for.
#

# In[15]:

res

# ## Step 4: Post-process, return result in classical format.
#

# If the procedure terminates correctly, then the prev_vector and iters values in our callback_dict dictionary should be equal to the solution vector and total number of function evaluations, respectively.  This is easy to verify:
#

# In[16]:

all(callback_dict["prev_vector"] == res.x)

# In[17]:

callback_dict["iters"] == res.nfev

# We can also now view the progress towards convergence as monitored by the cost history at each iteration:
#

# In[18]:

# Plot cost and classical distance history
fig, ax = plt.subplots(2, 1, figsize=(8, 8))
ax[0].plot(range(callback_dict["iters"]), callback_dict["cost_history"])
ax[0].set_xlabel("Iterations")
ax[0].set_ylabel("Cost")
ax[1].plot(range(callback_dict["iters"]), callback_dict["classical_distance_history"])
ax[1].set_xlabel("Iterations")
ax[1].set_ylabel("Classical Distance")
plt.tight_layout()
plt.show()

• Hey, welcome to QCSE! We prefer laser-focus questions, so that people having the same problem later on can then use the answer you got. In order for you to increase the odds of getting an answer, can you put your code directly in your post (there's a "code sample" button for that), and ask a single question per post? For instance, use this one to ask about changing x0, and post another question about the comparison to the classical method. Mar 17 at 20:42

If your main task is to run the TSP problem on a real hardware, then you don't have to change much in the original code, but keep in mind that in TSP number of variables increase by the square of number of nodes in the graph. For example if you are running a 5 node graph, the number of qubits required will be 25, and since the actual quantum backends are generally very busy, it might take some time for your code to actually run.

Regarding the code part follow the same tutorial till the

## Running it on Quantum Computer

part and then, if you have an account in IBMQ. Make sure to initiate it in the code file via:

from qiskit_ibm_runtime import QiskitRuntimeService, Sampler, Estimator, Options, Session

QiskitRuntimeService.save_account(channel="ibm_quantum", token="<Token ID here>",overwrite=True)
service = QiskitRuntimeService(channel='ibm_quantum')
backend = service.least_busy(operational=True, simulator=False)

print("The Backend is: ",backend.name)


and then for the sampler, do this:

options = Options()
options.transpilation.skip_transpilation = True
options.execution.shots = 100

session = Session(backend=backend)
sampler = Sampler(session=session,options=options)


and in the tutorial code, change this line:

vqe = SamplingVQE(sampler=Sampler(), ansatz=ry, optimizer=optimizer)


to this

vqe = SamplingVQE(sampler=sampler, ansatz=ry, optimizer=optimizer)


Also make sure to transpile your quantum circuit according to your backend using transpile function.

Run the rest as it is. You'll get the TSP result, but running on Quantum Hardware really takes a long time. If the queue is longer, it might even take days too. Just like this, it shows you what is the estimated time, it can be longer too.

• Thank you! I am now having a little trouble transpiling. I get the error AlgorithmError: 'The number of qubits of the ansatz does not match the operator, and the ansatz does not allow setting the number of qubits using num_qubits.' But I am using this code: from qiskit import transpile from qiskit.transpiler import PassManager transpiled_ry = transpile(ry, backend=backend, optimization_level=3) optimizer = SPSA(maxiter=300) vqe = SamplingVQE(sampler=Sampler(backend), ansatz=ansatz_isa, optimizer=optimizer) result = vqe.compute_minimum_eigenvalue(qubitOp) Mar 22 at 13:15
• Your circuit is saved in transpiled_ry, where does ansatz_isa come from? Mar 23 at 23:41
• Both I and my circuit did not make sense please let me try explain better. I tried transpiling: ry = TwoLocal(qubitOp.num_qubits, "ry", "cz", reps=5, entanglement="linear") transpiled_ry = transpile(ry, backend=backend, optimization_level=3) vqe = SamplingVQE(sampler=sampler, ansatz=transpiled_ry, optimizer=optimizer) result = vqe.compute_minimum_eigenvalue(qubitOp). I get the error AlgorithmError: 'The number of qubits of the ansatz does not match the operator, and the ansatz does not allow setting the number of qubits using num_qubits.' Mar 24 at 20:29
• So essentially, qubitOp which is my operator i.e. Hamiltonian for the TSP, and then transpiled_ry should be the transpiled ry waveform approximator quantum circuit. I thought that transpiling automatically sets up the qubits to work properly with the backend I am not sure why the numbers are mismatched. Specifically here, I have num_qubits = qubitOp.num_qubits ry = TwoLocal(qubitOp.num_qubits, "ry", "cz", reps=5, entanglement="linear") so num_qubits in ry and qubitOp should be the same, but qubitOp has 9 and transpiled_ry has 27 if I use num_qubits(). Mar 24 at 20:38
• Mar 25 at 4:10