DQI mistakes fixes (#1609)

A series of changes in the code to fix a mistake in the matrix

Author: daniela-angulo

demonstrations_v2/tutorial_dqi/demo.py

@@ -38,7 +38,7 @@
 where :math:`\mathbf{b_i}` is the :math:`i`-th row of matrix :math:`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
-:math:`v_i+\mathbf{b_i}\cdot\mathbf{x}` is always even, and that it is odd in the opposite case.
+:math:`v_i+\mathbf{b_i}\cdot\mathbf{x}` is always even, and 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 :math:`\mathbf{x}` from a uniform distribution.
@@ -104,7 +104,7 @@ def objective_function(x):
 #     :class: note
 #
 #     The `Hadamard transform <https://lucatrevisan.github.io/teaching/cs259q-12/lecture07b.pdf>`__ can be thought of as a version of the discrete Fourier transform. You have likely used it  
-#     several times when applying a Hadamard gate at the beginning of a circuit to prepare an equal superposition of all basis states. 
+#     several times when applying a layer of Hadamard gates at the beginning of a circuit to prepare an equal superposition of all basis states. 
 #     For :math:`n` qubits and a basis state :math:`|\mathbf{x}\rangle`, we can write the  
 #     transformation as
 # 
@@ -123,7 +123,7 @@ def objective_function(x):
 # where the coefficients :math:`w_k` are carefully chosen (see Section 8.1 of [#Jordan2024]_ for more about this result).
 # 
 # The DQI algorithm for solving the max-XORSAT problem involves three qubit registers: a weight, an
-# error, and a syndrome register, with dimensions :math:`\left \lceil \log_{2} \ell \right \rceil`,
+# error, and a syndrome register, with dimensions :math:`\left \lceil \log_{2} (\ell+1) \right \rceil`,
 # :math:`m`, and :math:`n`, respectively. The algorithm’s steps are outlined below and summarized in Figure 1.
 # 
 # 1. **Embed weight coefficients:** prepare the state :math:`\sum_{k=0}^{\ell} w_k|k\rangle` in the
@@ -161,10 +161,10 @@ def objective_function(x):
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # 
 # We are going to prepare the superposition :math:`\sum_{k=0}^{\ell} w_k|k\rangle`. As previously stated,
-# the weight register is made up of :math:`\left \lceil \log_{2} \ell \right \rceil` qubits, which means
+# the weight register is made up of :math:`\left \lceil \log_{2} (\ell+1) \right \rceil` qubits, which means
 # that the index :math:`k` is being binary encoded. The coefficients :math:`w_k` are chosen such that
 # they maximize the number of satisfied linear equations. These optimal weights are the components 
-# of the principal eigenvector of an :math:`(\ell+1)\times(\ell+1)` symmetric tridiagonal matrix [#Jordan2024]_:
+# of the principal eigenvector of an :math:`(\ell+1)\times(\ell+1)` symmetric tridiagonal matrix:
 # 
 # .. math::
 # 
@@ -183,8 +183,7 @@ def objective_function(x):
 # in this step.
 # 
 # In this demo, we will use a polynomial of degree :math:`2` for a reason that will become clear during
-# the decoding step. For now, you might be wondering if :math:`\left \lceil \log_{2} 2\right \rceil=1`
-# qubit will be enough to encode :math:`k=0,1,2`. Well, let’s examine what we obtain for the
+# the decoding step. For now, let’s examine what we obtain for the
 # coefficients vector :math:`\mathbf{w}`.
 # 
 
@@ -193,16 +192,16 @@ def objective_function(x):
 d = (p - 2 * r) / pnp.sqrt(r * (p - r))
 l = 2
 # Define registers
-num_weight_qubits = int(pnp.ceil(pnp.log2(l)))
+num_weight_qubits = int(pnp.ceil(pnp.log2(l+1)))
 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."""
+    """Calculates optimal weights: 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)
+    diag_sup = pnp.diag(pnp.sqrt(pnp.arange(1, l + 1) * (m - pnp.arange(1, l + 1) + 1)), 1)
     A = diag_main + diag_sup + pnp.transpose(diag_sup)
     _, eigenvectors = pnp.linalg.eigh(A)
     return eigenvectors[:, -1]
@@ -213,16 +212,14 @@ def w_k_optimal(m, l):
 
 
 ######################################################################
-# Since :math:`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
-# :math:`\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 :math:`k` values we are encoding are :math:`1` and :math:`2`.
+# We will need two qubits to encode :math:`w_k` for :math:`k=0,1,2`. The state is given by a weighted superposition
+# of the basis states :math:`|00\rangle, |01\rangle, |10\rangle`.
+# To prepare it, we will use ``qml.StatePrep`` in our ``embed_weights`` function shown below.
 #
 
 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)
+    """Prepares the weight register in superposition of given coefficients."""
+    qml.StatePrep(w_k, wires=weight_register, pad_with=0)
 
 
 ######################################################################
@@ -313,23 +310,28 @@ def format_state_vector(state_vector, tol: float = 1e-6):
 def SCS(m, k):
     """Implements the Split & Cycle shift unitary."""
 
+    # To address the correct set of wires
+    m_angle = m
+    m = m + num_weight_qubits - 1 
+
     # Two-qubit gate
     qml.CNOT(wires=[m - 1, m])
-    qml.CRY(2 * pnp.arccos(pnp.sqrt(1 / m)), wires=[m, m - 1])
+    qml.CRY(2 * pnp.arccos(pnp.sqrt(1 / m_angle)), 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.ctrl(qml.RY, (m, m - l + 1))(2 * pnp.arccos(pnp.sqrt(l / m_angle)), 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):
+    for wire_idx in range(m - k + num_weight_qubits, m + num_weight_qubits):
         qml.X(wires=wire_idx)
 
     # Applies the SCS unitaries
@@ -350,8 +352,8 @@ def weight_error_prep(m, n, l):
     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)
+    qml.ctrl(prepare_dicke_state, weight_register, control_values=(0,1))(m, k=1)
+    qml.ctrl(prepare_dicke_state, weight_register, control_values=(1,0))(m, k=2)
 
     return qml.state()
 
@@ -366,11 +368,11 @@ def weight_error_prep(m, n, l):
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # 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 :math:`m` with Hamming weight of
-# :math:`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 :math:`1`, as the qubit state was already :math:`|0\rangle`. From now on, we can 
-# choose to disregard the weight register.
+# We accomplished this by generating bit strings of length :math:`m` with Hamming weights of
+# :math:`1` and :math:`2` using the ``generate_bit_strings`` function. We then applied controlled bit flips to the
+# qubits in the :math:`1` state of the weight register. We did not need to perform any action for bit strings with
+# a Hamming weight of :math:`0`, as the qubit state was already :math:`|00\rangle`. From now on, we can 
+# disregard the weight register.
 # 
 
 from itertools import combinations
@@ -393,10 +395,19 @@ def generate_bit_strings(length, hamming_weight):
 
 
 def uncompute_weight(m, k):
-    """Uncomputes weight register when l=2"""
-    bit_strings = list(generate_bit_strings(m, k))
+    """Uncomputes weight register."""
+    bit_strings_dicke = list(generate_bit_strings(m, k))
+    binary_string_weight = bin(k)[2:].zfill(num_weight_qubits)
+
     for i in range(comb(m, k)):
-        qml.ctrl(qml.X, m_register, control_values=bit_strings[i])(0)
+        for j in range(len(binary_string_weight)):
+            bit = int(binary_string_weight[j])
+            if bit == 1:
+                qml.ctrl(
+                    qml.X, 
+                    m_register, 
+                    control_values=bit_strings_dicke[i]
+                )(weight_register[j])
 
 
 ######################################################################
@@ -405,24 +416,25 @@ def uncompute_weight(m, k):
 # 
 # To impart a phase :math:`(-1)^{\mathbf{v}\cdot\mathbf{y}}`, we perform a Pauli-Z on each qubit for
 # which :math:`v_i=1`. This is simply a conditional operation within a
-# ``for`` loop in the ``phase_Z`` function. Let’s now implement this step, together with the weight 
+# ``for`` loop in the ``phase_Z`` function. Let’s now implement this step, together with the weight register
 # uncomputation, in a function ``encode_v``, 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)
+        qml.cond(v[i],qml.Z)(wires=i + num_weight_qubits)
 
 
 @qml.qnode(dev)
 def encode_v(m, n, l):
     """Quantum circuit uncomputing weight register and encoding vector of constraints."""
 
     # Load the previous state
-    qml.StatePrep(raw_state_vector, wires=range(0, m + 1))
+    qml.StatePrep(raw_state_vector, wires=range(0, m + num_weight_qubits))
 
     # Uncompute weight register
+    uncompute_weight(m, k=1)
     uncompute_weight(m, k=2)
 
     # Impart phase
@@ -462,10 +474,10 @@ def B_T_multiplication(B_T, n_register):
 
 @qml.qnode(dev)
 def syndrome_prep(m, n, l):
-    """Quantum circuit preparing syndrome register"""
+    """Quantum circuit preparing syndrome register."""
 
     # Load the previous state
-    qml.StatePrep(raw_state_vector, wires=range(0, m + 1))
+    qml.StatePrep(raw_state_vector, wires=range(0, m + num_weight_qubits))
 
     # Compute s = B^T y into the syndrome register
     B_T_multiplication(B_T, n_register)
@@ -551,7 +563,7 @@ def decoding(m, n, l):
     """Quantum circuit decoding and uncomputing error register"""
 
     # Load the previous state
-    qml.StatePrep(raw_state_vector, wires=range(0, m + n + 1))
+    qml.StatePrep(raw_state_vector, wires=range(0, m + n + num_weight_qubits))
 
     # Uncompute syndrome register using a Lookup table
     for syndrome, error in decoding_table:
@@ -569,7 +581,7 @@ def decoding(m, n, l):
 # Hadamard transform and sample
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # 
-# After the previous step, we obtained the Hadamard transform of the state we are looking for. The final
+# After the previous step, we obtained the Hadamard transform of the state we were looking for. The final
 # step is to apply the Hadamard transform to this state to obtain
 # :math:`|P(f)\rangle=\sum_{\mathbf{x}} P(f(\mathbf{x}))|\mathbf{x}\rangle`. 
 # Let's write a ``DQI`` quantum function containing all the steps of the algorithm previously described. 
@@ -578,16 +590,17 @@ def decoding(m, n, l):
 @partial(qml.set_shots, shots=n_samples)
 @qml.qnode(dev)
 def DQI(m, n, l):
-    """Quantum circuit implementing DQI algorithm to solve max-XORSAT."""
+    """Quantum circuit implementing the DQI algorithm to solve max-XORSAT."""
 
     # Prepare weight register
-    qml.Hadamard(wires=0)
+    embed_weights(w_k,weight_register)
 
     # 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)
+    qml.ctrl(prepare_dicke_state, weight_register, control_values=(0,1))(m, k=1)
+    qml.ctrl(prepare_dicke_state, weight_register, control_values=(1,0))(m, k=2)
 
     # Uncompute weight register
+    uncompute_weight(m, k=1)
     uncompute_weight(m, k=2)
 
     # Impart phase
@@ -656,21 +669,27 @@ def DQI(m, n, l):
 # -----------
 #
 # .. [#Jordan2024]
+#
 #    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.
 #
 # .. [#Bartschi2019]
+#
 #    Andreas Bärtschi, and Stephan Eidenbenz.
 #    "Deterministic Preparation of Dicke States.",
 #    `<https://arxiv.org/abs/1904.07358>`__, 2019.
 #
 # .. [#Patamawisut2025]
+#
 #    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.
 #
-# .. [#Classiq2025] 
+# .. [#Classiq2025]
+# 
 #    Classiq. 
 #    “Decoded Quantum Interferometry Algorithm.”
 #    `<https://docs.classiq.io/latest/explore/algorithms/dqi/dqi_max_xorsat/>`__, 2025.
+#
+

demonstrations_v2/tutorial_dqi/metadata.json

@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2025-09-19T10:00:00+00:00",
-    "dateOfLastModification": "2025-09-22T15:48:14+00:00",
+    "dateOfLastModification": "2025-11-26T15:48:14+00:00",
     "categories": [
         "Algorithms"
     ],