Qubit-Efficient Variational Quantum Algorithms for Image Segmentation

🧠 Introduction

Classical algorithms are already incredibly effective at solving many optimization problems — especially when the number of variables is in the low thousands. However, when scaling up to problems involving tens or hundreds of thousands of binary variables, classical solvers begin to falter due to the combinatorial explosion in the search space. Quantum computing offers a compelling alternative — not by brute force, but by exploring exponentially large solution spaces in parallel using quantum superposition and entanglement. One popular quantum algorithm for such problems is the Quantum Approximate Optimization Algorithm (QAOA) ^[3], a textbook example of a variational quantum algorithm (VQA). However, QAOA requires one qubit per binary variable, which becomes prohibitively resource-intensive for large-scale problems. For instance, solving a QUBO (Quadratic Unconstrained Binary Optimization) problem with 10,000 variables using QAOA would require 10,000 logical qubits — a scale far beyond current quantum hardware capabilities.

🔍 Why Qubit-Efficient Methods?

This demo explores qubit-efficient alternatives to QAOA that are not only suitable for today's NISQ (Noisy Intermediate-Scale Quantum) devices, but also generalizable across combinatorial optimization problems.

We specifically demonstrate these methods in the context of unsupervised image segmentation, by:

Formulating segmentation as a min-cut problem over a graph derived from the image.
Reformulating the min-cut as a QUBO problem.
Solving the QUBO using qubit-efficient VQAs, namely ^[1]:
- Parametric Gate Encoding (PGE)
- Ancilla Basis Encoding (ABE)
- Adaptive Cost Encoding (ACE)

These methods only require a logarithmic number of qubits with respect to the problem size (i.e., number of binary variables). This makes them scalable and implementable even on near-term quantum hardware. Moreover, the exact graph-based image segmentation technique have been explored also using quantum annealers [2]. While this demo walks through an image segmentation example, the encoding schemes presented here can be applied to any binary combinatorial optimization problem.

Architecture for segmenting an image by transforming it into a graph and solving the corresponding minimum cut as a QUBO problem using variational quantum circuits ^[1].

Clone the Repository: Begin by cloning the repository to your local machine. Open a terminal/cmd/powershell and run

git clone https://github.com/supreethmv/Pennylane-ImageSegmentation.git
cd Pennylane-ImageSegmentation

Install Dependencies: Use the provided requirements.txt to install all necessary libraries. The process may take a few minutes depending on your internet speed and the packages already available in your system.

pip install -r requirements.txt

Open a jupyter notebook and execute the below python code.

🧱 Import dependencies

We begin by importing the standard Python libraries for graph construction and visualization, and PennyLane modules for quantum circuit construction and simulation. This setup ensures a seamless interface for hybrid quantum-classical processing.

import numpy as np
import networkx as nx
import pennylane as qml
from scipy.optimize import minimize, differential_evolution
import matplotlib.pyplot as plt

from math import ceil, log2
import time

🖼️ Generate an Example Image

To keep things intuitive, we generate a toy grayscale image. This synthetic image will be segmented using quantum algorithms.

We use a simple 4×4 grid for clarity. The pixel intensities will act as a proxy for similarity when we later construct a graph. Each pixel becomes a vertex, and edges represent similarity (e.g., inverse of intensity difference).

height,width = 4,4
image = np.array([
       [0.97,  0.92, 0.05, 0.92],
       [0.93,  0.10, 0.12, 0.93],
       [0.94,  0.15, 0.88, 0.94],
       [0.98,  0.15, 0.10, 0.92]
       ])


plt.imshow(image, interpolation='nearest', cmap=plt.cm.gray, vmin=0, vmax=1)
plt.xticks([])
plt.yticks([])
plt.show()

🔗 Representing Image as a Graph

Each image is modeled as an undirected weighted grid graph:

Each pixel becomes a node.
An edge exists between neighboring pixels.
The weight reflects similarity (e.g., inverse absolute difference in grayscale values).

This transforms the segmentation task into a minimum cut problem, where the goal is to partition the graph to minimize the sum of cut edge weights.

def image_to_grid_graph(gray_img, sigma=0.5):
  # Convert image to grayscale
  h, w = gray_img.shape
  # Initialize graph nodes and edges
  nodes = np.zeros((h*w, 1))
  edges = []
  nx_elist = []
  # Compute node potentials and edge weights
  min_weight = 1
  max_weight = 0
  for i in range(h*w):
    x, y = i//w, i%w
    nodes[i] = gray_img[x,y]
    if x > 0:
      j = (x-1)*w + y
      weight = 1-np.exp(-((gray_img[x,y] - gray_img[x-1,y])**2) / (2 * sigma**2))
      edges.append((i, j, weight))
      nx_elist.append(((x,y),(x-1,y),np.round(weight,2)))
      if min_weight>weight:min_weight=weight
      if max_weight<weight:max_weight=weight
    if y > 0:
      j = x*w + y-1
      weight = 1-np.exp(-((gray_img[x,y] - gray_img[x,y-1])**2) / (2 * sigma**2))
      edges.append((i, j, weight))
      nx_elist.append(((x,y),(x,y-1),np.round(weight,2)))
      if min_weight>weight:min_weight=weight
      if max_weight<weight:max_weight=weight
  a=-1
  b=1
  if max_weight-min_weight:
    normalized_edges = [(node1,node2,-1*np.round(((b-a)*((edge_weight-min_weight)/(max_weight-min_weight)))+a,2)) for node1,node2,edge_weight in edges]
    normalized_nx_elist = [(node1,node2,-1*np.round(((b-a)*((edge_weight-min_weight)/(max_weight-min_weight)))+a,2)) for node1,node2,edge_weight in nx_elist]
  else:
    normalized_edges = [(node1,node2,-1*np.round(edge_weight,2)) for node1,node2,edge_weight in edges]
    normalized_nx_elist = [(node1,node2,-1*np.round(edge_weight,2)) for node1,node2,edge_weight in nx_elist]
  return nodes, edges, nx_elist, normalized_edges, normalized_nx_elist

pixel_values, elist, nx_elist, normalized_elist, normalized_nx_elist = image_to_grid_graph(image)

pixel_values, normalized_nx_elist

(array([[0.97],
        [0.92],
        [0.05],
        [0.92],
        [0.93],
        [0.1 ],
        [0.12],
        [0.93],
        [0.94],
        [0.15],
        [0.88],
        [0.94],
        [0.98],
        [0.15],
        [0.1 ],
        [0.92]]),
 [((0, 1), (0, 0), np.float64(1.0)),
  ((0, 2), (0, 1), np.float64(-1.0)),
  ((0, 3), (0, 2), np.float64(-1.0)),
  ((1, 0), (0, 0), np.float64(1.0)),
  ((1, 1), (0, 1), np.float64(-0.9)),
  ((1, 1), (1, 0), np.float64(-0.92)),
  ((1, 2), (0, 2), np.float64(0.97)),
  ((1, 2), (1, 1), np.float64(1.0)),
  ((1, 3), (0, 3), np.float64(1.0)),
  ((1, 3), (1, 2), np.float64(-0.87)),
  ((2, 0), (1, 0), np.float64(1.0)),
  ((2, 1), (1, 1), np.float64(1.0)),
  ((2, 1), (2, 0), np.float64(-0.82)),
  ((2, 2), (1, 2), np.float64(-0.77)),
  ((2, 2), (2, 1), np.float64(-0.69)),
  ((2, 3), (1, 3), np.float64(1.0)),
  ((2, 3), (2, 2), np.float64(0.97)),
  ((3, 0), (2, 0), np.float64(1.0)),
  ((3, 1), (2, 1), np.float64(1.0)),
  ((3, 1), (3, 0), np.float64(-0.92)),
  ((3, 2), (2, 2), np.float64(-0.8)),
  ((3, 2), (3, 1), np.float64(1.0)),
  ((3, 3), (2, 3), np.float64(1.0)),
  ((3, 3), (3, 2), np.float64(-0.9))])

G = nx.grid_2d_graph(image.shape[0], image.shape[1])
G.add_weighted_edges_from(normalized_nx_elist)

def draw(G):
    plt.figure(figsize=(6,6))
    default_axes = plt.axes(frameon=True)
    pos = {(x,y):(y,-x) for x,y in G.nodes()}
    nx.draw_networkx(G, pos)
    nodes = nx.draw_networkx_nodes(G, pos, node_color=1-pixel_values,
                                  node_size=3000,
                                  cmap=plt.cm.Greys, vmin=0, vmax=1)
    nodes.set_edgecolor('k')
    edge_labels = nx.get_edge_attributes(G, "weight")
    nx.draw_networkx_edge_labels(G,
                                pos=pos,
                                edge_labels=edge_labels,
                                font_size=10,
                                font_color='tab:blue')
    plt.axis('off')

draw(G)
plt.show()

🔁 Parametric Gate Encoding (PGE)

In PGE, we directly encode the binary segmentation mask in the rotation parameters of a diagonal gate.

Only qubits are required, which is a huge gain over methods like QAOA that need $O(n)$.

Quantum Circuit:

Apply Hadamard gates to initialize a uniform superposition.
Apply a diagonal gate U(\vec{\theta}) whose entries encode binary values using a thresholding function $f(\theta_i)$.
Evaluate the cost function:

$$C(\vec{\theta}) = \langle \psi(\vec{\theta}) | L_G | \psi(\vec{\theta}) \rangle$$

where $L_G$ is the graph Laplacian.

The result is optimized using a classical optimizer, and $\theta_i$ are mapped to bits depending on whether $\theta_i < \pi$ or not.

n = G.number_of_nodes()
n_qubits = int(np.ceil(np.log2(n)))
n_qubits

def R(theta): 
    if abs(theta) > 2*np.pi or abs(theta) < 0:
        theta = abs(theta) - (np.floor(abs(theta)/(2*np.pi))*(2*np.pi))
    return 0 if 0 <= theta < np.pi else 1

# Define a PennyLane quantum device (using default.qubit simulator)
dev = qml.device("default.qubit", wires=n_qubits)

# Quantum circuit as a QNode
@qml.qnode(dev)
def circuit(params, H):
    # N = len(params)
    N = int(np.log2(len(params)))
    for i in range(N):
        qml.Hadamard(wires=i)

    diagonal_elements = [np.exp(1j * np.pi * R(param)) for param in params]
    # qml.DiagonalQubitUnitary(np.diag(diagonal_elements), wires=list(range(N)))
    qml.DiagonalQubitUnitary(diagonal_elements, wires=list(range(N)))

    return qml.expval(qml.Hermitian(H, wires=list(range(N))))

def cost_fn(params, H):
    global optimization_iteration_count
    optimization_iteration_count += 1
    cost = circuit(params, H)
    print(f'@ Iteration {optimization_iteration_count} Cost :', cost)
    return cost

def decode(optimal_params):
    return [R(param) for param in optimal_params]

def objective_value(x, w):
    return np.dot(x, np.dot(w, x))

def new_nisq_algo_solver(G, optimizer_method='Powell', initial_params_seed=123):
    global optimization_iteration_count
    optimization_iteration_count = 0
    n = G.number_of_nodes()
    w = nx.adjacency_matrix(G).todense()
    D_G = np.diag(list(dict(G.degree()).values()))
    A_G = w
    L_G = D_G - A_G
    n_padding = (2**int(np.ceil(np.log2(n)))-n)
    L_G = np.pad(L_G, [(0, n_padding), (0, n_padding) ], mode='constant')
    H = L_G

    np.random.seed(seed=initial_params_seed)
    initial_params = np.random.uniform(low=0.5, high=2*np.pi , size=(n+n_padding))

    result = minimize(
        fun=cost_fn,
        x0=initial_params,
        args=(H,),
        method=optimizer_method,
    )

    optimal_params, expectation_value = result.x, result.fun
    x = np.real(decode(optimal_params))
    x = x[:n]
    cut_value = objective_value(x, w)

    return x, expectation_value, cut_value

binary_solution, expectation_value, cut_value = new_nisq_algo_solver(G, optimizer_method = 'Powell')
binary_solution, expectation_value, cut_value

@ Iteration 1 Cost : 2.7137499999999983
@ Iteration 2 Cost : 2.7137499999999983
@ Iteration 3 Cost : 2.7137499999999983
@ Iteration 4 Cost : 2.7137499999999983
@ Iteration 5 Cost : 2.7137499999999983
.
.
.
@ Iteration 349 Cost : 0.18374999999999986
@ Iteration 350 Cost : 0.18374999999999986
@ Iteration 351 Cost : 0.18374999999999986
@ Iteration 352 Cost : 0.18374999999999986
@ Iteration 353 Cost : 0.18374999999999986

def decode_binary_string(x, height, width):
  mask = np.zeros([height, width])
  for index,segment in enumerate(x):
    mask[index//width,index%width]=segment
  return mask

mask = decode_binary_string(binary_solution, height, width)
plt.imshow(mask, interpolation='nearest', cmap=plt.cm.gray, vmin=0, vmax=1)
plt.xticks([])
plt.yticks([])
plt.show()

🧪 Ancilla Basis Encoding (ABE)

ABE uses $\log_2(n) + 1$ qubits:

$\log_2(n)$ register qubits represent binary states.
One ancilla qubit helps encode the binary decision.

Each computational basis state of the register corresponds to one binary variable. The ancilla qubit’s amplitude is used to decide the bit value ^[5]:

$$x_{v_i} = \begin{cases} 0 & \text{if } |a_i|^2 > |b_i|^2 \\ 1 & \text{otherwise} \end{cases}$$

The cost function is parameterized over expectation values of basis state projectors:

$$C(\vec{\theta}) = \sum_{i,j} Q_{ij} \frac{\langle \hat{P}_i^1 \rangle \langle \hat{P}_j^1 \rangle}{\langle \hat{P}i \rangle \langle \hat{P}j \rangle} (1 - \delta{ij}) + \sum_i Q{ii} \frac{\langle \hat{P}_i^1 \rangle}{\langle \hat{P}_i \rangle}$$

Example of a quantum circuit of ABE/ACE for solving QUBO of size 4 using 3 qubits (2 register + 1 ancilla), whereas in this demo we are solving a QUBO of size 16, so we are simulating 5 qubits. Basis state amplitudes are used to decode the binary solution (Fig. 3 in [1])

For the complete implementations of the next two methods ABE and ACE can be found in my personal blogging site available at:

https://www.supreethmv.com/blog/posts/qubit-efficient-encoding/

📄 References

[1] S. M. Venkatesh, A. Macaluso, M. Nuske, M. Klusch, and A. Dengel, "Qubit-Efficient Variational Quantum Algorithms for Image Segmentation," 2024 IEEE International Conference on Quantum Computing and Engineering (QCE), Montreal, QC, Canada, 2024, pp. 450–456. DOI: 10.1109/QCE60285.2024.00059, arXiv:2405.14405

[2] S. M. Venkatesh, A. Macaluso, M. Nuske, M. Klusch and A. Dengel, "Q-Seg: Quantum Annealing-Based Unsupervised Image Segmentation," in IEEE Computer Graphics and Applications, vol. 44, no. 5, pp. 27-39, Sept.-Oct. 2024, doi: 10.1109/MCG.2024.3455012, arXiv:2311.12912.

[3] Farhi, E., Goldstone, J., & Gutmann, S. (2014). A quantum approximate optimization algorithm. arXiv preprint arXiv:1411.4028.

[4] Rančić, M. J. (2023). Noisy intermediate-scale quantum computing algorithm for solving an n-vertex MaxCut problem with log (n) qubits. Physical Review Research, 5(1), L012021. DOI: 10.1103/PhysRevResearch.5.L012021

[5] Tan, B., Lemonde, M. A., Thanasilp, S., Tangpanitanon, J., & Angelakis, D. G. (2021). Qubit-efficient encoding schemes for binary optimisation problems. Quantum, 5, 454. DOI: 10.22331/q-2021-05-04-454

Qubit-Efficient Encoding Techniques for Solving QUBO Problems