- Demos/
- Algorithms/
Decoded Quantum Interferometry
Decoded Quantum Interferometry
Published: August 25, 2025. Last updated: August 25, 2025.
Showing quantum advantage is a long-sought goal by the quantum computing community. Clearly demonstrating that a quantum algorithm can solve a practical computational task with an advantage over classical algorithms remains an active area of research. Recently, Jordan et al. introduced an algorithm called Decoded quantum interferometry (DQI) to tackle combinatorial optimization problems 1. For a specific task of approximating optimal polynomial fits, it showed superior performance over existing classical counterparts. DQI takes a new approach; unlike QAOA and adiabatic optimization, it doesn’t use Hamiltonians to connect the optimization problem to quantum mechanics. Instead, it leverages quantum interference and transforms the optimization problem into a decoding one.
In this demo, we will implement DQI in PennyLane to solve the max-XORSAT problem. While this problem is simple and has not yet demonstrated quantum advantage, it clearly illustrates the operational principles and underlying intuition of DQI. We will begin by introducing the problem and the core principles of the algorithm, followed by a description of each algorithmic step and its implementation.
The max-XORSAT problem
Finding the best solution(s) from a large, finite set of possible candidates is known as combinatorial optimization. A well-known problem from this field is the traveling salesman problem. Mathematically, this can be phrased as trying to maximize an objective function whose domain is large and discrete.
The max-XORSAT problem is a simple example of this, where we are given an \(m\times n\) matrix \(B\) (with \(m>n\)) and a vector \(\mathbf{v}\) of length \(m\), and are required to find the \(n\)-bit string \(\mathbf{x}\) that satifies the maximum number of constraints imposed by the \(m\) linear mod-2 equations, over \(\mathbb{F}_2\), \(B\mathbf{x}=\mathbf{v}\). Here, \(\mathbf{v}\) is a vector of dimension \(m\) specified by the problem. Since this system of equations is over \(\mathbb{F}_2\), the matrix and vectors only contain zeros and ones.
The objective function we are aiming to maximize is
where \(\mathbf{b_i}\) is the \(i\)-th row of matrix \(B\). You can verify that this function represents the number of satisfied equations minus the number of unsatisfied ones by considering that when the equation is satisfied, the exponent \(v_i+\mathbf{b_i}\cdot\mathbf{x}\) is always even, and that it is odd in the opposite case.
Let’s define the conditions for our specific max-XORSAT problem and visualize the objective function in a histogram when randomly sampling bit strings \(\mathbf{x}\) from a uniform distribution. Later, we’ll use samples generated by the DQI algorithm and compare that objective function to this initial plot to see how well DQI performs.
import pennylane as qml
import math
from pennylane import numpy as pnp
import matplotlib.pyplot as plt
import itertools
from pprint import pprint
plt.style.use("pennylane.drawer.plot")
# Define parameters of max-XORSAT problem Bx=(mod2)v
B = pnp.array([[1, 0, 0, 0], [1, 1, 0, 0], [0, 1, 1, 0], [0, 0, 1, 1], [0, 0, 0, 1]])
v = pnp.array([1, 0, 1, 0, 1])
m, n = B.shape
B_T = B.T
n_samples = 500
def objective_function(x):
"""Calculates objective function."""
f = 0
for i in range(m):
f += (-1) ** (v[i] + pnp.dot(B[i], x))
return f
# Random sampling
samples = pnp.random.randint(0, 2, size=(n_samples, n))
f_x_array_random = [objective_function(sample) for sample in samples]
plt.hist(f_x_array_random, bins=30, density=True)
plt.xlabel(r"$f(x)$")
plt.ylabel("density")
plt.show()
The DQI Algorithm
Given the objective function presented earlier, a first approach to maximize it might be to prepare the state \(\sum_{\mathbf{x}} f(\mathbf{x})|\mathbf{x}\rangle\). This would increase the probability of sampling \(\mathbf{x}\) strings with high values of \(f(\mathbf{x})\). While this would work, DQI proposes a much more effective strategy to bias the sampling distribution and greatly enhance the probability of success. It proposes to encode a polynomial \(P(f(\mathbf{x}))\) of the function values in the amplitudes:
where \(P\) is of some degree \(\ell\). The challenge then becomes: how do we prepare such a state? DQI provides a concrete recipe.
The objective function \(f(\mathbf{x})\) has a very sparse Hadamard spectrum. This means that there are only \(m\) non-zero “frequency” components \(\mathbf{b_i}\) out of \(2^n\) possible in \(f(\mathbf{x})\). Consequently, preparing the state \(\sum_{\mathbf{x}} f(\mathbf{x})|\mathbf{x}\rangle\) can be seen as a simple task by creating an appropriate superposition of the \(m\) amplitudes and apply a Hadamard transform. The same principle holds for preparing \(|P(f)\rangle\), even though it is not as simple, we will still first prepare the Hadamard transform of the state, taking advantage of its sparse spectrum, and then transform back. The Hadamard transform of \(P(f(\mathbf{x}))\) is:
where the coefficients \(w_k\) are carefully chosen.
The DQI algorithm for solving the max-XORSAT problem involves three qubit registers: a weight, an error, and a syndrome register, with dimensions \(\left \lceil \log_{2} \ell \right \rceil\), \(m\), and \(n\), respectively. The algorithm’s steps are outlined below:
Embed weight coefficients: prepare the state \(\sum_{k=0}^{\ell} w_k|k\rangle\) in the weight register to encode the degree \(\ell\) of the polynomial.
Prepare Dicke states: generate Dicke states on the error register, conditioned on the value \(k\),
\[\begin{split}\sum_{k=0}^{\ell} w_k|k\rangle \frac{1}{\sqrt{\binom{m}{k}}}\sum_{\substack{\mathbf{y}\\|\mathbf{y}|=k}} |\mathbf{y}\rangle.\end{split}\]
Uncompute the weight register.
Encode the vector of constraints: encode \(\mathbf{v}\) by imparting a phase \((-1)^{\mathbf{v}\cdot\mathbf{y}}\),
\[\begin{split}\sum_{k=0}^{\ell} w_k \frac{1}{\sqrt{\binom{m}{k}}}\sum_{\substack{\mathbf{y}\\|\mathbf{y}|=k}} (-1)^{\mathbf{v}\cdot\mathbf{y}} |\mathbf{y}\rangle.\end{split}\]Compute syndrome: reversibly compute \(B^T \mathbf{y}\) into the syndrome register,
\[\begin{split}\sum_{k=0}^{\ell} w_k \frac{1}{\sqrt{\binom{m}{k}}}\sum_{\substack{\mathbf{y}\\|\mathbf{y}|=k}} (-1)^{\mathbf{v}\cdot\mathbf{y}} |\mathbf{y}\rangle|B^T \mathbf{y}\rangle.\end{split}\]Decode and uncompute: use the computed value of \(B^T \mathbf{y}\) to find \(\mathbf{y}\) and uncompute the error register. This results in the Hadamard transform of the desired state.
Hadamard transform and sample: apply the Hadamard transform to obtain \(|P(f)\rangle\) and sample from it to get solutions.
PennyLane implementation of DQI
Amplitude encoding: weight coefficients
We are going to prepare the superposition \(\sum_{k=0}^{\ell} w_k|k\rangle\). As previously stated, the weight register is made up of \(\left \lceil \log_{2} \ell \right \rceil\) qubits, which means that the index \(k\) is being binary encoded. The coefficients \(w_k\) are chosen such that they maximize the number of satisfied linear equations. The paper specifies that these optimal weights are the components of the principal eigenvector of an \((\ell+1)\times(\ell+1)\) symmetric tridiagonal matrix:
with \(a_k=\sqrt{k(m-k+1)}\) and \(d=\frac{p-2r}{\sqrt{r(p-r)}}\). Here, \(p\) is the number of elements of the finite field where our problem lives (in this case, \(p=2\)), and \(r\) is the number of inputs that will yield \(f=+1\) (for this problem, \(r=1).\) In this demo, we will use a polynomial of degree \(2\) for a reason that will become clear during the decoding step. For now, you might be wondering if \(\left \lceil \log_{2} 2\right \rceil=1\) qubit will be enough to encode \(k=0,1,2\). Well, let’s examine what we obtain for the coefficients vector \(\mathbf{w}\).
p = 2
r = 1
d = (p - 2 * r) / pnp.sqrt(r * (p - r))
l = 2
# Define registers
num_weight_qubits = int(pnp.ceil(pnp.log2(l)))
weight_register = range(num_weight_qubits)
m_register = range(num_weight_qubits, m + num_weight_qubits)
n_register = range(m + num_weight_qubits, n + m + num_weight_qubits)
def w_k_optimal(m, l):
"""Calculates optimal weights for superposition: principal vector of tridiagonal matrix."""
diag_main = pnp.diag(pnp.arange(l + 1) * d)
diag_sup = pnp.diag(pnp.sqrt(pnp.arange(l) * (m - pnp.arange(l) + 1)), 1)
A = diag_main + diag_sup + pnp.transpose(diag_sup)
_, eigenvectors = pnp.linalg.eigh(A)
return eigenvectors[:, -1]
w_k = w_k_optimal(m, l)
print("the optimal values for w are", w_k)
the optimal values for w are [0. 0.70710678 0.70710678]
Since \(w_0=0\), a single qubit is sufficient to encode the remaining non-zero coefficients. The
explicit form of this state will be the uniform superposition
\(\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)\), which can be readily prepared by a Hadamard gate. However,
to make our code more versatile, we will use qml.StatePrep in our embed_weights function shown below.
For subsequent steps.Let’s just keep in mind that the \(k\) values we are encoding are \(1\) and \(2\).
def embed_weights(w_k, weight_register):
"""Prepares the weight register to be in superposition given coefficients."""
qml.StatePrep(w_k[1:], wires=weight_register, pad_with=0)
Prepare Dicke states with k excitations
For this step, we need a conditional operation that prepares Dicke states \(|D^{m}_{k}\rangle\)—where the superscript is the number of qubits and the subscript is the number of excitations—for each index \(k\) in the weight register. For our particular example, we will prepare Dicke states with one and two excitations.
Before implementing the conditional operation, let’s briefly review a method for preparing such states as presented in 2. This deterministic algorithm provides a unitary \(U_{m,k}\) to generate a Dicke state \(|D_k^m\rangle\) from \(|0\rangle^{\otimes m-k}|1\rangle^{\otimes k}\) as input. The method is based on the fact that Dicke states can be expressed in an inductive form as
Given that decomposition, one can generate the two Dicke states on the right-hand side of the equation by applying a smaller unitary \(U_{m-1,k}\) provided a suitable superposition of the inputs \(|0\rangle^{\otimes m-q}|1\rangle^{\otimes q-1}\) and \(|0\rangle^{\otimes m-q}|1\rangle^{\otimes q}\) is prepared first. Following the paper’s convention, we define a Split and Cyclic shift unitary \(\mathrm{SCS}_{m,k}\) to prepare such suitable superposition:
This inductive decomposition implies that \(U_{m,k}\) can be implemented by applying \(\mathrm{SCS}_{m,k}\) followed by the smaller unitary \(U_{m-1,k}\). This process continues; we keep decomposing the unitaries into a split and cycle operation followed by a smaller unitary until we reach the base case \(U_{1,1}=\mathrm{Id}\). Ultimately, \(U_{m,k}\) is composed of a series of subsequently smaller \(\mathrm{SCS}_{m,k}\) operations, as summarized by the following equation:
The explicit form of \(\mathrm{SCS}_{m,k}\), which is implemented in code by the function
SCS, is not discussed here. However, it’s worth noting that it is composed of a two-qubit gate
followed by \(k-1\) three-qubit gates (further details can be found in 2).
Let’s now implement this algorithm in the prepare_dicke_state function and then use it twice in
a controlled way via qml.ctrl() in our main function DQI where all the parts of the DQI
algorithm are going to be placed. To verify that this step was implemented correctly, we will output
the quantum state of the weight and error registers printed in a nice form using the auxiliary format_state_vector function.
def format_state_vector(state_vector, tol: float = 1e-6):
"""Formats a state vector as a dictionary of bit-strings and amplitudes."""
num_qubits = int(pnp.log2(len(state_vector)))
state_dict = {}
for i, amplitude in enumerate(state_vector):
if pnp.abs(amplitude) > tol:
# Format the index 'i' as a binary string with leading zeros
bit_string = format(i, f'0{num_qubits}b')
state_dict[bit_string] = complex(amplitude)
return state_dict
def SCS(m, k):
"""Implements the Split & Cycle shift unitary."""
# Two-qubit gate
qml.CNOT(wires=[m - 1, m])
qml.CRY(2 * pnp.arccos(pnp.sqrt(1 / m)), wires=[m, m - 1])
qml.CNOT(wires=[m - 1, m])
# k-1 three-qubit gates
for l in range(2, k + 1):
qml.CNOT(wires=[m - l, m])
qml.ctrl(qml.RY, (m, m - l + 1))(2 * pnp.arccos(pnp.sqrt(l / m)), wires=m - l)
qml.CNOT(wires=[m - l, m])
def prepare_dicke_state(m, k):
"""Prepares a Dicke state with m qubits and k excitations in a inductive form."""
# Prepares input state
for wire_idx in range(m - k + 1, m + 1):
qml.X(wires=wire_idx)
# Applies the SCS unitaries
for i in reversed(range(k + 1, m + 1)):
SCS(i, k)
for i in reversed(range(2, k + 1)):
SCS(i, i - 1)
dev = qml.device("default.qubit")
@qml.qnode(dev)
def DQI(m, n, l):
"""Quantum circuit implementing DQI algorithm to solve max-XORSAT."""
# Prepare weight register
embed_weights(w_k, weight_register)
# Prepare Dicke states conditioned on k values
qml.ctrl(prepare_dicke_state, (0,), control_values=0)(m, k=1)
qml.ctrl(prepare_dicke_state, (0,), control_values=1)(m, k=2)
return qml.state()
raw_state_vector = DQI(m, n, l)
formatted_state = format_state_vector(raw_state_vector)
pprint(formatted_state)
{'000001': (0.316227766016838+0j),
'000010': (0.31622776601683783+0j),
'000100': (0.3162277660168379+0j),
'001000': (0.31622776601683794+0j),
'010000': (0.3162277660168379+0j),
'100011': (0.2236067977499789+0j),
'100101': (0.22360679774997896+0j),
'100110': (0.22360679774997896+0j),
'101001': (0.223606797749979+0j),
'101010': (0.223606797749979+0j),
'101100': (0.22360679774997896+0j),
'110001': (0.22360679774997896+0j),
'110010': (0.22360679774997896+0j),
'110100': (0.2236067977499789+0j),
'111000': (0.22360679774997896+0j)}
Uncompute the weight register
After preparing the Dicke states, we uncompute and discard the state of the weight register. In
general, this is a straightforward process, as the Hamming weights encoded are known.
We accomplished this by generating bit strings of length \(m\) with Hamming weight of
\(2\) using the generate_bit_strings function. We then applied a controlled bit flip to the
weight register for these specific cases. We did not need to perform any action for bit strings with
a Hamming weight of \(1\), as the qubit state was already \(|0\rangle\). From now on, we can
choose to discard it and not include it in our outputs.
def generate_bit_strings(length, hamming_weight):
"""Generates all bit strings of a given length and Hamming weight."""
one_positions = itertools.combinations(range(length), hamming_weight)
results = []
for positions in one_positions:
# Create a new list of zeros
bit_string = pnp.zeros(length, dtype=int)
# Place 1s at the specified positions
bit_string[pnp.array(positions)] = 1
results.append(bit_string.tolist())
return results
def uncompute_weight(m, k):
"""Uncomputes weight register when l=2"""
bit_strings = list(generate_bit_strings(m, k))
for i in range(math.comb(m, k)):
qml.ctrl(qml.X, m_register, control_values=bit_strings[i])(0)
Encode constraints vector
To impart a phase \((-1)^{\mathbf{v}\cdot\mathbf{y}}\), we perform a Pauli-Z on each qubit for
which \(v_i=1\), as instructed by the DQI paper. This is simply a conditional operation within a
for loop in the phase_Z function. Let’s now include this step in the DQI function and
output the resulting quantum state.
def phase_Z(v):
"""Imparts a phase (-1)^{vy}."""
for i in range(len(v)):
qml.cond(v[i],qml.Z)(wires=i + 1)
dev = qml.device("default.qubit")
@qml.qnode(dev)
def DQI(m, n, l):
"""Quantum circuit implementing DQI algorithm to solve max-XORSAT."""
# Prepare weight register
qml.Hadamard(wires=0)
# Prepare Dicke states conditioned on k values
qml.ctrl(prepare_dicke_state, (0), control_values=0)(m, 1)
qml.ctrl(prepare_dicke_state, (0), control_values=1)(m, 2)
# Uncompute weight register
uncompute_weight(m, k=2)
# Impart phase
phase_Z(v)
return qml.state()
raw_state_vector = DQI(m, n, l)
formatted_state = format_state_vector(raw_state_vector)
pprint(formatted_state)
{'000001': (-0.316227766016838+0j),
'000010': (0.31622776601683783+0j),
'000011': (-0.2236067977499789+0j),
'000100': (-0.3162277660168379+0j),
'000101': (0.22360679774997896-0j),
'000110': (-0.22360679774997896+0j),
'001000': (0.31622776601683794+0j),
'001001': (-0.223606797749979+0j),
'001010': (0.223606797749979+0j),
'001100': (-0.22360679774997896+0j),
'010000': (-0.3162277660168379+0j),
'010001': (0.22360679774997896-0j),
'010010': (-0.22360679774997896+0j),
'010100': (0.2236067977499789-0j),
'011000': (-0.22360679774997896+0j)}
Encode matrix B in the syndrome register
We are almost finished, hang in there! Now, we need to compute \(B^T \mathbf{y}\) in the syndrome register. While it may not be immediately obvious how to implement this as a unitary operation, the binary nature of the matrix and vector allows for a smooth translation. This operation can be realized using CNOT gates, with controls on the error register and targets on the syndrome register 3. Specifically, a CNOT is applied for every entry \(B_{ij}^T=1\), controlled on the \(j\)-th qubit of the error register and with the \(i\)-th qubit of the syndrome register as the target.
def B_T_multiplication(B_T, n_register):
"""Computes B^T y into the syndrome register."""
for row_index, row in enumerate(B_T):
for col_index, element in enumerate(row):
if element == 1:
qml.CNOT(wires=[m_register[col_index], n_register[row_index]])
dev = qml.device("default.qubit", shots=n_samples)
@qml.qnode(dev)
def DQI(m, n, l):
"""Quantum circuit implementing DQI algorithm to solve max-XORSAT."""
# Prepare weight register
qml.Hadamard(wires=0)
# Prepare Dicke states conditioned on k values
qml.ctrl(prepare_dicke_state, (0), control_values=0)(m, 1)
qml.ctrl(prepare_dicke_state, (0), control_values=1)(m, 2)
# Uncompute weight register
uncompute_weight(m, k=2)
# Impart phase
phase_Z(v)
# Compute s = B^T y into the syndrome register
B_T_multiplication(B_T, n_register)
return qml.counts(wires=range(1, m + n + 1))
pprint(DQI(m, n, l))
{'000010001': 46,
'000100011': 51,
'000110010': 29,
'001000110': 47,
'001010111': 24,
'001100101': 21,
'010001100': 50,
'010011101': 29,
'010101111': 26,
'011001010': 31,
'100001000': 52,
'100011001': 15,
'100101011': 28,
'101001110': 28,
'110000100': 23}
Decoding!
This step is the main challenge of the algorithm: uncomputing the error register \(\mathbf{y}\) using the information from the syndrome register \(\mathbf{s}=B^T \mathbf{y}\) in an efficient way. This would be an easy task if \(B\) were always a square matrix; however, since it is not, we need to solve an underdetermined linear system of equations \(\mathbf{s}=B^T \mathbf{y}\) subject to the constraint \(|\mathbf{y}|\leq \ell\) given by the known Hamming weights of \(\mathbf{y}\). The problem we have just described is precisely the syndrome decoding problem, where \(B^T\) is the parity-check matrix, \(\mathbf{s}\) is the syndrome, and \(\mathbf{y}\) is the error.
Polynomial degree \(\ell\)
The kernel of \(B^T\) defines an error-correcting code. The distance \(d\) of this code determines the number of errors it can correct, given by \(\left \lfloor (d-1)/2 \right \rfloor\). This condition ensures that the decoding problem has a unique solution. For this demo, we choose \(\ell=2\) such that it is less than half the distance of the code \(d=5\) and this condition is met. For a detailed discussion of the restrictions on \(\ell\), please refer to the original paper 1.
To keep things simple, we will use a straightforward approach for decoding by building a Lookup Table
(LUT) in which we compute the syndrome for each possible error using the classical function
syndrome_LUT 4. Then, for each syndrome in the syndrome register, the corresponding error is
uncomputed in the error register using controlled bit-flip operations. We will now integrate this
into ourDQI function and see in the output how the syndrome register is uncomputed.
def syndrome_LUT(parity_check_matrix_T):
"""Generates the Lookup table given a parity-check matrix."""
num_data_qubits = parity_check_matrix_T.shape[1]
syndrome_dict = {}
for i in range(2**num_data_qubits):
error_bitstring = format(i, f"0{num_data_qubits}b")
error_vector = pnp.array([int(bit) for bit in error_bitstring])
syndrome_vector = pnp.dot(parity_check_matrix_T, error_vector) % 2
syndrome_bitstring = "".join(map(str, syndrome_vector))
if syndrome_bitstring not in syndrome_dict:
syndrome_dict[syndrome_bitstring] = error_bitstring
else:
existing_error = syndrome_dict[syndrome_bitstring]
if pnp.sum(error_vector) < pnp.sum(pnp.array([int(bit) for bit in existing_error])):
syndrome_dict[syndrome_bitstring] = error_bitstring
# Convert the dictionary to the desired list-of-lists format
lookup_matrix = []
for syndrome_str, error_str in syndrome_dict.items():
syndrome_list = [int(bit) for bit in syndrome_str]
error_list = [int(bit) for bit in error_str]
lookup_matrix.append([syndrome_list, error_list])
# Sort the matrix for a consistent output (e.g., by syndrome)
lookup_matrix.sort(key=lambda x: int("".join(map(str, x[0])), 2))
return lookup_matrix
# Generate the lookup table
decoding_table = syndrome_LUT(B_T)
dev = qml.device("default.qubit", shots=n_samples)
@qml.qnode(dev)
def DQI(m, n, l):
"""Quantum circuit implementing DQI algorithm to solve max-XORSAT."""
# Prepare weight register
qml.Hadamard(wires=0)
# Prepare Dicke states conditioned on k values
qml.ctrl(prepare_dicke_state, (0), control_values=0)(m, 1)
qml.ctrl(prepare_dicke_state, (0), control_values=1)(m, 2)
# Uncompute weight register
uncompute_weight(m, k=2)
# Impart phase
phase_Z(v)
# Compute s = B^T y into the syndrome register
B_T_multiplication(B_T, n_register)
# Uncompute syndrome register using a Lookup table
for syndrome, error in decoding_table:
for i in range(len(error)):
if error[i] == 1:
qml.ctrl(qml.X, n_register, control_values=syndrome)(m_register[i])
return qml.counts(wires=range(1, m + n + 1))
pprint(DQI(m, n, l))
{'000000001': 52,
'000000010': 25,
'000000011': 57,
'000000100': 21,
'000000101': 20,
'000000110': 58,
'000000111': 21,
'000001000': 44,
'000001001': 21,
'000001010': 28,
'000001011': 28,
'000001100': 50,
'000001101': 26,
'000001110': 25,
'000001111': 24}
Hadamard transform and sample
After the previous step, we obtained the Hadamard transform of the state were looking for. The final step is to apply the Hadamard transform to this state to obtain \(|P(f)\rangle=\sum_{\mathbf{x}} P(f(\mathbf{x}))|\mathbf{x}\rangle\). Finally, we will collect samples from this state, calculate their objective values, and build a histogram to compare with the random sampling done at first.
dev = qml.device("default.qubit", shots=n_samples)
@qml.qnode(dev)
def DQI(m, n, l):
"""Quantum circuit implementing DQI algorithm to solve max-XORSAT."""
# Prepare weight register
qml.Hadamard(wires=0)
# Prepare Dicke states conditioned on k values
qml.ctrl(prepare_dicke_state, (0), control_values=0)(m, 1)
qml.ctrl(prepare_dicke_state, (0), control_values=1)(m, 2)
# Uncompute weight register
uncompute_weight(m, k=2)
# Impart phase
phase_Z(v)
# Compute s = B^T y into the syndrome register
B_T_multiplication(B_T, n_register)
# Uncompute syndrome register using a Lookup table
for row in decoding_table:
syndrome, error = row
for i in range(len(error)):
if error[i] == 1:
qml.ctrl(qml.X, n_register, control_values=syndrome)(m_register[i])
# Apply Hadamard transform
for wire in n_register:
qml.Hadamard(wire)
return qml.counts(wires=n_register)
results = DQI(m, n, l)
all_bit_lists = []
for string, samples in results.items():
integer_list = [int(bit) for bit in string]
for _ in range(samples):
all_bit_lists.append(integer_list)
DQI_array = pnp.array(all_bit_lists) # DQI sample array
f_x_array_DQI = [objective_function(DQI_array[i]) for i in range(n_samples)]
pprint(results)
plt.hist(f_x_array_random, bins=30, alpha=0.7, density=True, label="random sampling")
plt.hist(f_x_array_DQI, bins=30, alpha=0.5, density=True, label="DQI sampling")
plt.xlabel(r"$f(x)$")
plt.ylabel("density")
plt.legend(prop={"size": 10})
plt.show()
{'0000': 23,
'0001': 20,
'0010': 21,
'0011': 51,
'0100': 29,
'0101': 17,
'0110': 10,
'0111': 23,
'1000': 24,
'1001': 21,
'1010': 19,
'1011': 51,
'1100': 60,
'1101': 58,
'1110': 20,
'1111': 53}
The histogram shows a significantly higher probability of obtaining a bit string with a high objective value when sampling from the DQI algorithm compared to uniform sampling. It is also worth noting that the samples produced by DQI are unbiased, meaning that the probability is uniform for all solutions within a given objective value.
Conclusion
In this tutorial, we implemented the Decoded Quantum Interferometry algorithm to solve the max-XORSAT problem. We’ve shown how DQI recasts a combinatorial optimization task as a decoding problem by leveraging quantum interference and encoding a polynomial of the objective function in the amplitudes of a quantum state. While we used a simplified problem here, the principles we’ve explored are applicable to more complex optimization tasks. This algorithm marks the beginning of an exploration into a different approach for achieving quantum speedups.
References
- 1(1,2)
Stephen P. Jordan, Noah Shutty, Mary Wootters, Adam Zalcman, Alexander Schmidhuber, Robbie King, Sergei V. Isakov, Tanuj Khattar, and Ryan Babbush. “Optimization by Decoded Quantum Interferometry.”, https://arxiv.org/abs/2408.08292, 2024.
- 2(1,2)
Andreas Bärtschi, and Stephan Eidenbenz. “Deterministic Preparation of Dicke States.”, https://arxiv.org/abs/1904.07358, 2019.
- 3
Natchapol Patamawisut, Naphan Benchasattabuse, Michal Hajdušek, and Rodney Van Meter. “Quantum Circuit Design for Decoded Quantum Interferometry.” https://arxiv.org/abs/2504.18334, 2025.
- 4
Classiq. “Decoded Quantum Interferometry Algorithm.” https://docs.classiq.io/latest/explore/algorithms/dqi/dqi_max_xorsat/, 2025.
Total running time of the script: (0 minutes 1.042 seconds)
