Solving the Quadratic Assignment Problem : A Classically Hard Benchmark for Quantum Optimization

How quantum algorithms like VQE, CVaR-VQE handle constrained optimization problems too tough for classical solvers.

Introduction: The Hardest of the Hard

The Quadratic Assignment Problem (QAP) is widely regarded as one of the most difficult problems in combinatorial optimization, often labeled “hardest of the hard” due to its resistance to classical methods even at moderate sizes.

First introduced by Koopmans and Beckmann in 1957, QAP models the task of assigning n facilities to n locations in a way that minimizes the total cost, which depends on both the flow between facilities and the distance between locations.

Despite decades of algorithmic development, solving QAP instances beyond 20–30 nodes remains a formidable challenge for even the most advanced solvers. Its structure combines:

• Dense quadratic cost interactions (via flow × distance terms), and
• Hard permutation constraints (each facility must go to exactly one location).

This makes QAP a natural test-bed to push the boundaries of quantum optimization algorithms. In this article, we explore how the QAP can be encoded and solved using emerging quantum methods, treating it not just as a toy problem, but as a true benchmark for scaling quantum advantage.

Why QAP Is Exceptionally Hard

QAP isn’t just NP-hard, it sits at the top of the difficulty hierarchy in practical optimization. Some key reasons:

• Quadratic explosion: Unlike linear assignment or knapsack problems, the cost function in QAP involves all pairwise interactions. This makes naive enumeration scale as O(n!) with dense cost structure.
• Permutation space: Valid solutions must satisfy strict bijective mappings (assignments), forming a highly constrained feasible set.
• Non-convexity: The objective is typically non-convex, non-separable, and riddled with local minima, making local search ineffective.
• Solver bottlenecks: Classical solvers like Gurobi, CPLEX, or SCIP often time out or return high optimality gaps even for instances with n≈25.

In bench-marking challenges (like QAPLIB), QAP consistently emerges as one of the last frontiers for optimization solvers, a perfect storm of computational hardness.

Why QAP Is Ideal for Quantum Benchmarks

The properties that make QAP classically difficult are the same features that test the limits of quantum optimization algorithms:

Unlike toy problems (e.g. Max-Cut on sparse graphs), QAP pushes both the qubit count and the circuit fidelity required to represent realistic combinatorial instances. Solving QAP provides meaningful evidence of progress.

Where is QAP Used?

The Quadratic Assignment Problem may sound abstract, but it models a wide range of real-world problems that require assigning entities to locations while minimizing interaction costs. It appears in domains such as:

• Facility Layout Planning: Deciding where to place departments in a factory to minimize transportation of materials based on flow and distance.
• Hospital Design: Assigning rooms or units to locations to reduce patient transfer times and improve efficiency.
• Keyboard Design: Determining the optimal layout of keys to minimize finger movement based on character transition frequency.
• Parallel Task Scheduling: Mapping computational tasks to processors where communication cost depends on task dependencies and processor distances.
• Data Center Optimization: Assigning virtual machines to physical servers to reduce communication latency and energy usage.

These applications share a common trait: interactions between pairs matter, and assignments affect overall cost non-linearly. That’s exactly what QAP captures — making it not only theoretically hard, but also practically valuable.

Formulating QAP as QUBO

The QAP can be encoded into a Quadratic Unconstrained Binary Optimization (QUBO) model by defining binary variables:

The objective becomes:

With additional constraints ensuring each facility is assigned to exactly one location, and each location gets exactly one facility, either through penalty terms or constraint-preserving mixers (in QAOA).

How we formulate and solve the QAP?

We are going to make use of QAPLIB bench-marking instances, provided here: https://coral.ise.lehigh.edu/data-sets/qaplib/qaplib-problem-instances-and-solutions/

A sample problem instance looks like:

12

    0    90    10    23    43     0     0     0     0     0    0     0
   90     0     0     0     0    88     0     0     0     0    0     0
   10     0     0     0     0     0    26    16     0     0    0     0
   23     0     0     0     0     0     0     0     0     0    0     0
   43     0     0     0     0     0     0     0     0     0    0     0
    0    88     0     0     0     0     0     0     1     0    0     0
    0     0    26     0     0     0     0     0     0     0    0     0
    0     0    16     0     0     0     0     0     0    96    0     0
    0     0     0     0     0     1     0     0     0     0   29     0
    0     0     0     0     0     0     0    96     0     0    0    37
    0     0     0     0     0     0     0     0    29     0    0     0
    0     0     0     0     0     0     0     0     0    37    0     0

    0    36    54    26    59    72     9    34    79    17   46    95
   36     0    73    35    90    58    30    78    35    44   79    36
   54    73     0    21    10    97    58    66    69    61   54    63
   26    35    21     0    93    12    46    40    37    48   68    85
   59    90    10    93     0    64     5    29    76    16    5    76
   72    58    97    12    64     0    96    55    38    54    0    34
    9    30    58    46     5    96     0    83    35    11   56    37
   34    78    66    40    29    55    83     0    44    12   15    80
   79    35    69    37    76    38    35    44     0    64   39    33
   17    44    61    48    16    54    11    12    64     0   70    86
   46    79    54    68     5     0    56    15    39    70    0    18
   95    36    63    85    76    34    37    80    33    86   18     0

To load this data, load basic packages :

import numpy as np
from docplex.mp.model import Model

and since the data files are .dat files, we need to extract data from them, using the following function:

def parse_qap_dat_file(file_path):
    """
    Parses a .dat file for the Quadratic Assignment Problem (QAP).

    Parameters:
    - file_path: str, path to the .dat file.

    Returns:
    - n: int, the size of the problem (number of facilities/locations).
    - A: 2D numpy array, flow matrix.
    - B: 2D numpy array, distance matrix.
    """
    with open(file_path, 'r') as file:
        lines = file.readlines()

    # Read the size of the problem
    n = int(lines[0].strip())

    # Read the flow matrix (A)
    A = []
    current_line = 1
    for i in range(n):
        row = list(map(int, lines[current_line].strip().split()))
        while len(row) < n:  # Handle cases where rows are split across multiple lines
            current_line += 1
            row.extend(list(map(int, lines[current_line].strip().split())))
        A.append(row[:n])
        current_line += 1

    # Read the distance matrix (B)
    B = []
    for i in range(n):
        row = list(map(int, lines[current_line].strip().split()))
        while len(row) < n:  # Handle cases where rows are split across multiple lines
            current_line += 1
            row.extend(list(map(int, lines[current_line].strip().split())))
        B.append(row[:n])
        current_line += 1

    return n, np.array(A), np.array(B)

def create_qap_model(n, A, B):
    """
    Creates a CPLEX model for the Quadratic Assignment Problem (QAP).

    Parameters:
    - n: int, the size of the problem (number of facilities/locations).
    - A: 2D numpy array, flow matrix.
    - B: 2D numpy array, distance matrix.

    Returns:
    - model: CPLEX model.
    - x: 2D list of CPLEX binary variables representing the assignment.
    """
    # Create a CPLEX model
    model = Model(name="Quadratic Assignment Problem")

    # Decision variables: x[i][j] = 1 if facility i is assigned to location j
    x = [[model.binary_var(name=f"x_{i}_{j}") for j in range(n)] for i in range(n)]

    # Objective: Minimize the total cost
    model.minimize(
        model.sum(A[i, k] * B[j, l] * x[i][j] * x[k][l]
                  for i in range(n) for j in range(n) for k in range(n) for l in range(n))
    )

    # Constraints: Each facility is assigned to exactly one location
    for i in range(n):
        model.add_constraint(model.sum(x[i][j] for j in range(n)) == 1, f"facility_assignment_{i}")

    # Constraints: Each location is assigned exactly one facility
    for j in range(n):
        model.add_constraint(model.sum(x[i][j] for i in range(n)) == 1, f"location_assignment_{j}")

    return model, x

and you can select the path:

file_path = "qapdata/chr12a.dat"

and then you can extract individual parameters:

n, A, B = parse_qap_dat_file(file_path)

# Create the QAP model
model, x = create_qap_model(n, A, B)
print(model.export_to_string())

you’ll get something like

\ This file has been generated by DOcplex
\ ENCODING=ISO-8859-1
\Problem name: Quadratic Assignment Problem

Minimize
 obj: [ 12960 x_0_0*x_1_1 + 19440 x_0_0*x_1_2 + 9360 x_0_0*x_1_3
      + 21240 x_0_0*x_1_4 + 25920 x_0_0*x_1_5 + 3240 x_0_0*x_1_6
      + 12240 x_0_0*x_1_7 + 28440 x_0_0*x_1_8 + 6120 x_0_0*x_1_9
      + 16560 x_0_0*x_1_10 + 34200 x_0_0*x_1_11 + 1440 x_0_0*x_2_1
      + 2160 x_0_0*x_2_2 + 1040 x_0_0*x_2_3 + 2360 x_0_0*x_2_4
      + 2880 x_0_0*x_2_5 + 360 x_0_0*x_2_6 + 1360 x_0_0*x_2_7 + 3160 x_0_0*x_2_8
      + 680 x_0_0*x_2_9 + 1840 x_0_0*x_2_10 + 3800 x_0_0*x_2_11
      + 3312 x_0_0*x_3_1 + 4968 x_0_0*x_3_2 + 2392 x_0_0*x_3_3
      + 5428 x_0_0*x_3_4 + 6624 x_0_0*x_3_5 + 828 x_0_0*x_3_6 + 3128 x_0_0*x_3_7
      + 7268 x_0_0*x_3_8 + 1564 x_0_0*x_3_9 + 4232 x_0_0*x_3_10
      + 8740 x_0_0*x_3_11 + 6192 x_0_0*x_4_1 + 9288 x_0_0*x_4_2
      + 4472 x_0_0*x_4_3 + 10148 x_0_0*x_4_4 + 12384 x_0_0*x_4_5
      + 1548 x_0_0*x_4_6 + 5848 x_0_0*x_4_7 + 13588 x_0_0*x_4_8
      + 2924 x_0_0*x_4_9 + 7912 x_0_0*x_4_10 + 16340 x_0_0*x_4_11
      + 12960 x_0_1*x_1_0 + 26280 x_0_1*x_1_2 + 12600 x_0_1*x_1_3
      + 32400 x_0_1*x_1_4 + 20880 x_0_1*x_1_5 + 10800 x_0_1*x_1_6
      + 28080 x_0_1*x_1_7 + 12600 x_0_1*x_1_8 + 15840 x_0_1*x_1_9
      + 28440 x_0_1*x_1_10 + 12960 x_0_1*x_1_11 + 1440 x_0_1*x_2_0
      + 2920 x_0_1*x_2_2 + 1400 x_0_1*x_2_3 + 3600 x_0_1*x_2_4
      + 2320 x_0_1*x_2_5 + 1200 x_0_1*x_2_6 + 3120 x_0_1*x_2_7
      + 1400 x_0_1*x_2_8 + 1760 x_0_1*x_2_9 + 3160 x_0_1*x_2_10
      + 1440 x_0_1*x_2_11 + 3312 x_0_1*x_3_0 + 6716 x_0_1*x_3_2
      + 3220 x_0_1*x_3_3 + 8280 x_0_1*x_3_4 + 5336 x_0_1*x_3_5
      + 2760 x_0_1*x_3_6 + 7176 x_0_1*x_3_7 + 3220 x_0_1*x_3_8
      + 4048 x_0_1*x_3_9 + 7268 x_0_1*x_3_10 + 3312 x_0_1*x_3_11
      + 6192 x_0_1*x_4_0 + 12556 x_0_1*x_4_2 + 6020 x_0_1*x_4_3
......
....

which is the LP of the problem.

And, to solve it with classical methods, like CPLEX, you can:

solution = model.solve()

if solution:
    print("Optimal value (total cost):", solution.objective_value)
    assignment = [(i, j) for i in range(n) for j in range(n) if x[i][j].solution_value > 0.5]
    print("Optimal assignment:", assignment)
else:
    print("No solution found.")

Optimal value (total cost): 9552.0
Optimal assignment: [(0, 6), (1, 4), (2, 11), (3, 1), (4, 0), (5, 2), (6, 8), (7, 10), (8, 9), (9, 5), (10, 7), (11, 3)]

This is a relatively small problem, hence it will be easier to solve it via CPLEX.

Converting QAP to QUBO

To solve the Quadratic Assignment Problem (QAP) on a quantum computer, we must first translate it into a Quadratic Unconstrained Binary Optimization (QUBO) formulation. This step is crucial because quantum algorithms like QAOA and quantum annealing operate on QUBO or Ising models.

Step 1: Define Binary Variables

We introduce binary variables:

For n facilities and n locations, we need n² binary variables. A valid solution must assign each facility to one unique location, and each location must receive exactly one facility, forming a permutation matrix.

🧠 Step 2: Encode the Objective

The classical QAP objective is:

This is naturally quadratic in terms of binary variables — making it directly suitable as the QUBO objective function.

🧷 Step 3: Add Constraint Penalties

Because QUBO doesn’t allow hard constraints, we introduce penalty terms to ensure each facility is assigned to one location and vice versa.

One-hot row constraint (each facility to one location):

One-hot column constraint (each location gets one facility):

These are converted to soft penalties by adding quadratic terms to the objective:

Final QUBO Form

The final QUBO objective becomes:

Doing it via the code, we load the basic Qiskit packages:

from qiskit_optimization.translators import from_docplex_mp
from qiskit_optimization.converters import QuadraticProgramToQubo

and making the quadratic program:

qp = from_docplex_mp(model)
print(qp.export_as_lp_string())

\ This file has been generated by DOcplex
\ ENCODING=ISO-8859-1
\Problem name: Quadratic Assignment Problem

Minimize
 obj: [ 12960 x_0_0*x_1_1 + 19440 x_0_0*x_1_2 + 9360 x_0_0*x_1_3
      + 21240 x_0_0*x_1_4 + 25920 x_0_0*x_1_5 + 3240 x_0_0*x_1_6
      + 12240 x_0_0*x_1_7 + 28440 x_0_0*x_1_8 + 6120 x_0_0*x_1_9
      + 16560 x_0_0*x_1_10 + 34200 x_0_0*x_1_11 + 1440 x_0_0*x_2_1
      + 2160 x_0_0*x_2_2 + 1040 x_0_0*x_2_3 + 2360 x_0_0*x_2_4
      + 2880 x_0_0*x_2_5 + 360 x_0_0*x_2_6 + 1360 x_0_0*x_2_7 + 3160 x_0_0*x_2_8
      + 680 x_0_0*x_2_9 + 1840 x_0_0*x_2_10 + 3800 x_0_0*x_2_11
      + 3312 x_0_0*x_3_1 + 4968 x_0_0*x_3_2 + 2392 x_0_0*x_3_3
      + 5428 x_0_0*x_3_4 + 6624 x_0_0*x_3_5 + 828 x_0_0*x_3_6 + 3128 x_0_0*x_3_7
      + 7268 x_0_0*x_3_8 + 1564 x_0_0*x_3_9 + 4232 x_0_0*x_3_10
      + 8740 x_0_0*x_3_11 + 6192 x_0_0*x_4_1 + 9288 x_0_0*x_4_2
      + 4472 x_0_0*x_4_3 + 10148 x_0_0*x_4_4 + 12384 x_0_0*x_4_5
      + 1548 x_0_0*x_4_6 + 5848 x_0_0*x_4_7 + 13588 x_0_0*x_4_8
      + 2924 x_0_0*x_4_9 + 7912 x_0_0*x_4_10 + 16340 x_0_0*x_4_11
      + 12960 x_0_1*x_1_0 + 26280 x_0_1*x_1_2 + 12600 x_0_1*x_1_3
      + 32400 x_0_1*x_1_4 + 20880 x_0_1*x_1_5 + 10800 x_0_1*x_1_6
      + 28080 x_0_1*x_1_7 + 12600 x_0_1*x_1_8 + 15840 x_0_1*x_1_9
      + 28440 x_0_1*x_1_10 + 12960 x_0_1*x_1_11 + 1440 x_0_1*x_2_0
      + 2920 x_0_1*x_2_2 + 1400 x_0_1*x_2_3 + 3600 x_0_1*x_2_4
      + 2320 x_0_1*x_2_5 + 1200 x_0_1*x_2_6 + 3120 x_0_1*x_2_7....

and then making it to QUBO:

converter = QuadraticProgramToQubo()
qubo = converter.convert(qp)

print(qubo.export_as_lp_string())

\ This file has been generated by DOcplex
\ ENCODING=ISO-8859-1
\Problem name: Quadratic Assignment Problem

Minimize
 obj: - 23823940 x_0_0 - 23823940 x_0_1 - 23823940 x_0_2 - 23823940 x_0_3
      - 23823940 x_0_4 - 23823940 x_0_5 - 23823940 x_0_6 - 23823940 x_0_7
      - 23823940 x_0_8 - 23823940 x_0_9 - 23823940 x_0_10 - 23823940 x_0_11
      - 23823940 x_1_0 - 23823940 x_1_1 - 23823940 x_1_2 - 23823940 x_1_3
      - 23823940 x_1_4 - 23823940 x_1_5 - 23823940 x_1_6 - 23823940 x_1_7
      - 23823940 x_1_8 - 23823940 x_1_9 - 23823940 x_1_10 - 23823940 x_1_11
      - 23823940 x_2_0 - 23823940 x_2_1 - 23823940 x_2_2 - 23823940 x_2_3
      - 23823940 x_2_4 - 23823940 x_2_5 - 23823940 x_2_6 - 23823940 x_2_7
      - 23823940 x_2_8 - 23823940 x_2_9 - 23823940 x_2_10 - 23823940 x_2_11
      - 23823940 x_3_0 - 23823940 x_3_1 - 23823940 x_3_2 - 23823940 x_3_3
      - 23823940 x_3_4 - 23823940 x_3_5 - 23823940 x_3_6 - 23823940 x_3_7
      - 23823940 x_3_8 - 23823940 x_3_9 - 23823940 x_3_10 - 23823940 x_3_11
      - 23823940 x_4_0 - 23823940 x_4_1 - 23823940 x_4_2 - 23823940 x_4_3
      - 23823940 x_4_4 - 23823940 x_4_5 - 23823940 x_4_6 - 23823940 x_4_7
      - 23823940 x_4_8 - 23823940 x_4_9 - 23823940 x_4_10 - 23823940 x_4_11.....

...
..
.

We can see the number of variables using:

# number of variables
num_vars = qubo.get_num_vars()
print('Number of variables:', num_vars)

Number of variables: 144

and then following the same steps as we did in the previous articles:

• MIS — [https://medium.com/@_monitsharma/solving-the-maximum-independent-set-problem-a-classically-hard-benchmark-for-quantum-optimization-5d8e9c267052]
• MDKP — [https://medium.com/@_monitsharma/solving-the-multi-dimensional-knapsack-problem-a-classically-hard-benchmark-for-quantum-2610e955d114]

or trying any quantum algorithm after that, taken from :

• SMU Quantum Algorithms — [ https://github.com/SMU-Quantum/quantum-optimization-algorithms]

or any problem instance:

• SMU Quantum Benchmarks — [https://github.com/SMU-Quantum/quantum-optimization-benchmarks/tree/main]

Comparing Classical and Quantum Performance on QAP

To assess how quantum algorithms perform on the QAP, we compare them directly against classical baselines across multiple dimensions:

As you can see, most methods just fail to give a feasible solution, while the others who give a solution are really bad.

While the time taken by them was very high:

Try it Yourself

All of the code, instances, and quantum formulations used in this article are available in our open-source GitHub repository:

🔗 GitHub: Quantum Optimization Benchmarks

also in Medium-Articles

The QAP experiments live inside this folder: 📁 quantum-optimization-benchmarks /Quadratic_Assignment_Problem/

You’ll find:

• Ready-to-run QUBO instances
• Scripts for VQE and Pauli Correlation Encoding
• Jupyter notebooks for interactive exploration

Read the full paper : A Comparative Study of Quantum Optimization Techniques for Solving Combinatorial Optimization Benchmark Problems

Repo with all code and instances : SMU-Quantum

You can also ⭐️ star the repo to be notified when new instances, plots, and articles are released.

Stay tuned — and follow the project for updates:

📝 Medium: @_MonitSharma 💻 GitHub: MonitSharma 🐦 Twitter: @_MonitSharma

🔜 Coming Next: Quantum Optimization Faces Market Share Allocation

This article is the next in an n-part series investigating how quantum algorithms scale on classically hard benchmark problems in discrete optimization.

Our next focus: the Market Share Problem, a combinatorial challenge that models how to allocate products to retailers while minimizing deviation from a target market distribution.

We’ll explore:

• How to formulate the Market Share Problem as a QUBO, capturing the exact and approximate fulfillment of retailer demands
• Why this problem is hard: it reduces to the NP-complete subset sum problem when simplified, and generalizes to a binary feasibility problem over hyperplane intersections
• What makes the classical version brittle, especially in high-dimensional, tightly constrained cases with conflicting retailer targets
• How quantum approaches such as QAOA, CVaR-VQE perform in terms of feasibility, objective gap, and solution reliability

This is where quantum optimization must go beyond objective maximization, it must find structure in sparse constraints, tolerate infeasibility, and still get as close to the market targets as possible.

The Market Share Problem isn’t just about profit — it’s about precision under pressure.