From Hello World to MaxCut: Solving Graph Problems with Qiskit

From Basic Principles to Real-World Applications in Graph Optimization

Previously we covered quick introduction of qiskit. Qiskit team provides a good course to help advance from basics to real-world problems and patterns to solve them.

Getting started with Qiskit

Qiskit 101 & 102

Key concepts

Backend: Represents the quantum hardware (a real quantum computer) or a simulator. Common backends: IBM Quantum devices, simulators like Aer (a simulator in Qiskit).

Primitives: smallest processing instruction for a given abstraction level. For primitives to work as effectively as possible, they should ideally be supported by both Qiskit and the backend provider. The Qiskit module qiskit.primitives provides the required support on the Qiskit side, and providers like IBM’s Qiskit Runtime enable access to appropriate backends through native implementations of their own primitives.

Quantum circuit: A sequence of quantum operations applied to qubits.

Circuit library: collection of valuable, well-studied circuits and gates (e.g. phase estimation, Fourier transform, or variational circuit) that Qiskit users can easily plug into their quantum circuits. These include useful tools for a wide array of quantum computing applications — including quantum machine learning, quantum chemistry, circuits for benchmarking.

Observables: In quantum mechanics, observables correspond to physical properties that can be measured (expectation value of an operator, used as input to loss/objective function). When considering a system of spins, for example, you could be interested in measuring the system’s energy or obtaining information about the alignment of the spins, such as the magnetization or the correlations between spins. To measure an n-qubit observable O on a quantum computer, you must represent it as a sum of tensor products of Pauli operators {I, X, Y, Z} (how quantum states evolve is unitary). Prove Pauli operators are unitary with code below

import numpy as np
from qiskit import QuantumCircuit
from qiskit.quantum_info import Operator

# Function to check unitarity
def check_unitary(gate):
    # Get the operator corresponding to the gate
    op = Operator(gate)
    # Compute the adjoint (conjugate transpose)
    adjoint = op.adjoint()
    # Check if U * U† = I
    product = op @ adjoint
    identity = np.eye(op.dim[0])  # Identity matrix of the same dimension
    # Returns True if two arrays are element-wise equal within a tolerance
    return np.allclose(product, identity)
# Define single-qubit gates
def create_hadamard_gate():
    circ_h = QuantumCircuit(1)
    circ_h.h(0)
    return circ_h.to_gate()
def create_pauli_x_gate():
    circ_x = QuantumCircuit(1)
    circ_x.x(0)
    return circ_x.to_gate()
def create_pauli_y_gate():
    circ_y = QuantumCircuit(1)
    circ_y.y(0)
    return circ_y.to_gate()
def create_pauli_z_gate():
    circ_z = QuantumCircuit(1)
    circ_z.z(0)
    return circ_z.to_gate()
def create_s_gate():
    circ_s = QuantumCircuit(1)
    circ_s.s(0)
    return circ_s.to_gate()
def create_t_gate():
    circ_t = QuantumCircuit(1)
    circ_t.t(0)
    return circ_t.to_gate()
    
# Dictionary of gates
gates = {
    "Hadamard": create_hadamard_gate(),
    "Pauli X": create_pauli_x_gate(),
    "Pauli Y": create_pauli_y_gate(),
    "Pauli Z": create_pauli_z_gate(),
    "S Gate": create_s_gate(),
    "T Gate": create_t_gate(),
}
# Check unitarity for each gate
for gate_name, gate in gates.items():
    is_unitary = check_unitary(gate)
    print(f"{gate_name} is unitary: {is_unitary}")

Estimator: computes expectation values of observables with respect to states prepared by quantum circuits. Users generally specify a list of circuits, observables, and possibly some additional configuration, with which the program can efficiently evaluate expectation values and variances. Workflow: Initialize the account, Create a circuit and an observable and transformed to only use instructions supported by the QPU (referred to as instruction set architecture (ISA) circuits), Initialize Qiskit Runtime Estimator, Invoke the Estimator and get results

Sampler: sampling from quantum circuits’ outputs, calculates probabilities or quasi-probabilities of bitstrings from quantum circuits (raw outcomes (the bitstrings measured in the computational basis)). Workflow: Initialize the account, Create a circuit and transformed to only use instructions supported by the QPU (referred to as instruction set architecture (ISA) circuits), Initialize the Qiskit Runtime Sampler, Invoke the Sampler and get results

Qiskit 102

To run the notebook, you need GraphViz installed

brew install graphviz
pip3 install graphviz

Key concepts

Hamiltonian: a function used to solve a problem of optimal control for a dynamical system.

Hamiltonian simulation: a problem that demands algorithms which implement the evolution of a quantum state efficiently.

Hamiltonian (quantum mechanics): The Hamiltonian of a system represents the total energy of the system; that is, the sum of the kinetic and potential energies of all particles associated with the system. Here, the problem is mapped to a Hamiltonian whose ground state corresponds to the optimal solution (objective function).

Quantum-classical hybrid algorithms: considering the noisy intermediate-scale quantum (NISQ) era, they are leveraged to overcome current limitations. Quantum parallelism is part of the strategy, but the overall algorithm is not fully parallel in the way you might expect from purely quantum algorithms like Grover’s. Generally, quantum parallelism exists in the sense that the quantum state can encode many possible solutions simultaneously in a superposition of states. When measured, the quantum state collapses, and the result corresponds to one of these possible solutions, with the probability distribution determined by the quantum state. However, the process of searching for the optimal solution or performing optimization is still iterative, as the quantum part provides an estimate of the objective function, and the classical optimizer updates the quantum circuit parameters based on this feedback.

Ising Models: It serves as a natural interface between physics and computer science, particularly in mapping to quadratic unconstrained binary optimization (QUBO) models. It was a model for magnetism, A physical system can be completely described by giving an expression for its energy. An Ising model consists of “spins” which can take values of +1 or −1. The energy of an Ising model is written in terms of whether spins are the same or different, which can be written as products of these terms. The mathematic equation can be leveraged by combinatorial optimization problems.

Quantum Approximate Optimization Algorithm (QAOA): one of the most representative quantum-classical hybrid algorithms, is designed to solve combinatorial optimization problems by transforming the discrete optimization problem (e.g. maxcut, traveling salesman, knapsack problem, QUBO) into a classical optimization problem over continuous circuit parameters. In QAOA, the quantum circuit is parameterized by continuous variables (angles of quantum gates (gamma for cost Hamiltonian and beta for mixer Hamiltonian (often a simple operator like the sum of Pauli-X operators acting on each qubit, which encourages exploration of different states))), which are then optimized classically to find an approximate solution to the discrete combinatorial problem.

General workflow of QAOA

1. State Initialization: common starting state is a uniform superposition of all possible bitstrings (representing all possible solutions)
2. Ansatz (quantum circuit), apply cost Hamiltonian and mixer Hamiltonian evolution (repeat for p layers)
3. Parameter optimization optimizes the variational parameters gamma and beta using a classical optimizer by evaluating the cost function on the quantum computer for different parameter settings and updating the parameters iteratively to minimize the cost function
4. Repeat or solution extraction measures the final state of the QAOA circuit to obtain an approximate solution to the optimization problem. Depending on the problem, you may also run multiple rounds of QAOA to get a statistical distribution of solutions and choose the best one.

Using Maxcut graph problem, which partitions graph’s vertices into two complementary sets S and T, such that the number of edges between S and T is as large as possible. Objective function:

• Objective: Find an assignment of xi∈{+1,−1} for all vertices ii that maximizes C(x).
• Maximizing C(x): Ensures that as many edges as possible (or edges with maximum total weight) are “cut” by the partition, i.e., their endpoints are in different subsets.

Steps of modeling optimization problem

Steps of modeling optimization problem in qiskit, with quantum-classical method. It is somewhat similar to machine learning/neural network, 2.

1. Setup objective function (like loss function in ML) that represents the problem goal
2. Setup model architecture for the problem
3. Initialize model and iterate
4. Setup real objective function that leverages objective function and model to retrieve value from hardware
5. Setup optimizer

Setup objective function that represents the problem goal

Generally in quantum computing, classical optimization problems are represented as quantum Hamiltonians whose ground states (lowest energy states) correspond to optimal solutions of the original problem. When we measure an observable, the possible outcomes are the eigenvalues of its corresponding operator. Since xi∈{+1,−1}, Pauli Z operator Z^i has eigenvalues ±1, xi​→Z^i​ (we can use Z^i to replace xi​ in objective function). Also since in quantum mechanics, systems naturally evolve towards states of minimum energy, we change maximize to minimize. Therefore, our Hamiltonian for maxcut is

which is 1/2 * (number of edges) — 1/2 * sum(Z^i​ Z^j), which are the terms.

def max_cut_hamiltonian(graph: rx.PyGraph) -> SparsePauliOp:
    # "" represents identity operator, [] indicates acting on no qubits
    terms = [["", [], 0.5 * graph.num_edges()]]
    for edge in graph.edge_list():
        # "ZZ": Pauli Z operators act on two qubits.
        # edge: A tuple (i, j) specifying the qubits (vertices) index operators act upon
        # -0.5: coefficient of the term
        terms.append(["ZZ", edge, -0.5])
    return SparsePauliOp.from_sparse_list(terms, num_qubits=graph.num_nodes())

Setup model architecture for the problem

We need to setup cost Hamiltonian and mixer Hamiltonian.

cost Hamiltonian: related to the above objective function, but needs to be expressed in quantum gate implementation. Since above objective function main term is ZZ (on pair of qubits (qubit i and qubit j) of an edge), here we use RZZGate. Angle gamma controls the strength of the interaction between qubits.

mixer Hamiltonian: Pauli X operator acting on qubit i. This mixer term encourages exploration of the solution space by rotating the qubits. Angle beta controls the strength of the mixing applied to each qubit (the mixer unitary facilitates transitions between different possible solutions (bitstrings)/changes probability amplitudes among the states, which enhances the amplitudes of good solutions (constructive interference) and suppress bad ones (destructive interference). It helps the quantum state to move from one potential solution to another and prevents the algorithm from getting stuck in local minima by allowing it to explore a broader region of the solution space).

def qaoa_circuit(graph: rx.PyGraph, p: int) -> QuantumCircuit:
    betas = [Parameter(f"b_{i}") for i in range(p)]
    gammas = [Parameter(f"g_{i}") for i in range(p)]
    # each qubit is mapped to one node in the graph
    qubits = QuantumRegister(graph.num_nodes())
    circuit = QuantumCircuit(qubits)
    circuit.h(qubits)
    for beta, gamma in zip(betas, gammas):
        for (i, j) in graph.edge_list():
            circuit.append(RZZGate(gamma), [qubits[i], qubits[j]])
        for q in qubits:
            circuit.append(RXGate(beta), [q])
    return circuit

# corresponding to 5 nodes in the graph
n_qubits = 5
probability = 0.8
# layer
p = 1
graph = rx.undirected_gnp_random_graph(num_nodes=n_qubits, probability=probability)
observable = max_cut_hamiltonian(graph)
circuit = qaoa_circuit(graph, p=p)

Initialize model and iterate

Remember we have two parameters in our circuit: gamma and beta. The following code initializes them.

from qiskit_ibm_runtime import EstimatorV2 as Estimator
import numpy as np
rng = np.random.default_rng()
params = rng.uniform(-np.pi, np.pi, size=circuit.num_parameters)

estimator = Estimator(mode=backend)
pub = (isa_circuit, isa_observable, params)
job = estimator.run([pub])

Setup real objective function

def f(x: np.ndarray, estimator: Estimator) -> float:
    pub = (isa_circuit, isa_observable, x)
    job = estimator.run([pub])
    result = job.result()
    pub_result = result[0]
    val = -float(pub_result.data.evs)
    print(f"Objective function value: {val}")
    return val

Setup optimizer

The sample uses scipy.optimize.minimize as the optimizer. Main arguments include: objective function to be minimized, initial parameters/guess, extra arguments passed to the objective function (e.g. estimator to compute expectation value from circuit).

from qiskit_ibm_runtime import Session
import scipy.optimize

with Session(backend=backend):
    estimator = Estimator(mode=backend)
    # max 5 iterations
    result = scipy.optimize.minimize(
        f, params, args=(estimator,), method="COBYLA", options=dict(maxiter=5)
    )

Solve utility-scale quantum optimization problems

Ansatze and Variational Forms

The above article talks about another way to create Hamiltonian and QAOAAnsatz for a graph. However, this way requires you to provide the mapping from graph edge to operator explicitly, e.g.

For following graph

from qiskit.circuit.library import QAOAAnsatz
from qiskit.quantum_info import SparsePauliOp

# Problem to Hamiltonian operator
hamiltonian = SparsePauliOp.from_list([("ZZII", 1), ("IZZI", 1), ("ZIIZ", 1), ("IZIZ", 1), ("IIZZ", 1)])
# QAOA ansatz circuit
ansatz = QAOAAnsatz(hamiltonian, reps=2)

Explain for SparsePauliOp arguments. Remember the 0th qubit is farthest right, so edge from node 2 to 3 maps to ZZII (from right to left, II means no operation on qubit 0 and 1, ZZ means operation on qubit 2 and 3). Similarly, IZZI refers to the edge between node 1 and 2. Since there is no edge between node 1 and 3, so there is no “ZIZI”.

To prove they are same, we can create a graph with rustnetworkx and use hamiltonian function above and compare

import rustworkx as rx
from rustworkx.visualization import mpl_draw

n = 4
graph = rx.PyGraph()
graph.add_nodes_from(range(n))
# The edge syntax is (start, end, weight)
edges = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0)]
graph.add_edges_from(edges)

mpl_draw(graph, pos=rx.shell_layout(graph), with_labels=True, edge_labels=str, node_color="#1192E8")

nx_observable = max_cut_hamiltonian(graph)
nx_observable

Result (since no edge between node 1 and 3, it adds “IIII”. Others are same between 2 approaches.

SparsePauliOp(['IIII', 'IIZZ', 'IZIZ', 'ZIIZ', 'IZZI', 'ZZII'],

SparsePauliOp(['ZZII', 'IZZI', 'ZIIZ', 'IZIZ', 'IIZZ'],

Bonus

CHSH Inequality

Use quantum computing to validate Bell Inequalities (CHSH inequality is a specific form of Bell’s inequality), which is used to test whether the correlations predicted by quantum mechanics can be explained by local hidden variables (local realism).

The tutorial is good to follow, with the following complementary explanations.

The observables are related to

S=E(AB)−E(Ab)+E(aB)+E(ab)

with following conventions:

two qubits representing two parties

• Alice: A and a
• A corresponds to Z (Pauli-Z on the first qubit).
• a corresponds to X (Pauli-X on the first qubit).
• Bob: B and b
• B corresponds to Z (Pauli-Z on the second qubit).
• b corresponds to X (Pauli-X on the second qubit).

Arguments of SparsePauliOp

(“ZX”, -1): meaning operating on two qubits (no need to pass qubit index explicitly), Pauli-Z on the first qubit, Pauli-X on the second qubit, coefficient is -1 (matching −E(Ab))

# <CHSH1> = <AB> - <Ab> + <aB> + <ab> -> <ZZ> - <ZX> + <XZ> + <XX>
observable1 = SparsePauliOp.from_list([("ZZ", 1), ("ZX", -1), ("XZ", 1), ("XX", 1)])

# <CHSH2> = <AB> + <Ab> - <aB> + <ab> -> <ZZ> + <ZX> - <XZ> + <XX>
observable2 = SparsePauliOp.from_list([("ZZ", 1), ("ZX", 1), ("XZ", -1), ("XX", 1)])

The theta parameter in the circuit is essential because it enables you to simulate the different measurement settings required to test the CHSH inequality. Quantum measurements are basis-dependent. The default measurement in a quantum circuit is in the computational (Z) basis. To measure in a different basis (e.g., at an angle θ), you need to rotate the qubit so that the new basis aligns with the computational basis. The ry(theta) gate achieves this rotation. In CHSH experiment, Alice and Bob (the two observers) independently choose between different measurement settings (angles) for their qubits. The ry(theta) gate models Alice’s choice of measurement setting. By varying θ, you can simulate different measurement angles that Alice might choose.

theta = Parameter("$\\theta$")

chsh_circuit = QuantumCircuit(2)
chsh_circuit.h(0)
chsh_circuit.cx(0, 1)
chsh_circuit.ry(theta, 0)
chsh_circuit.draw(output="mpl", idle_wires=False, style="iqp")