Mthree

qadence_protocols/mitigations/readout.py

@@ -6,10 +6,35 @@
 import numpy.typing as npt
 from numpy.linalg import inv, matrix_rank, pinv
 from qadence import QuantumModel
+from qadence.logger import get_logger
 from qadence.noise.protocols import Noise
-from qadence.types import ReadOutOptimization
 from scipy.linalg import norm
 from scipy.optimize import LinearConstraint, minimize
+from scipy.sparse import csr_matrix
+from scipy.sparse.linalg import gmres
+
+from qadence_protocols.types import ReadOutOptimization
+
+logger = get_logger(__name__)
+
+
+def normalized_subspace_kron(noise_matrices: npt.NDArrary, subspace: npt.NDArray) -> npt.NDArray:
+    n_qubits = len(noise_matrices)
+    conf_matrix = csr_matrix((2**n_qubits, 2**n_qubits))
+
+    for j in subspace:
+        for i in subspace:
+            bin_i = bin(i)[2:].zfill(n_qubits)
+            bin_j = bin(j)[2:].zfill(n_qubits)
+
+            # Manually computing the entries of tensor product for only the subspace
+            conf_matrix[i, j] = np.prod(
+                [noise_matrices[k][int(bin_i[k])][int(bin_j[k])] for k in range(n_qubits)]
+            )
+
+        conf_matrix[:, j] /= np.sum(conf_matrix[:, j])
+
+    return conf_matrix
 
 
 def tensor_rank_mult(qubit_ops: npt.NDArray, prob_vect: npt.NDArray) -> npt.NDArray:
@@ -26,7 +51,7 @@ def tensor_rank_mult(qubit_ops: npt.NDArray, prob_vect: npt.NDArray) -> npt.NDAr
 
     # Contract each tensor index (qubit) with the inverse of the single-qubit
     for i in range(N):
-        prob_vect_t = np.tensordot(qubit_ops[i], prob_vect_t, axes=(1, i))
+        prob_vect_t = np.tensordot(qubit_ops[N - 1 - i], prob_vect_t, axes=(1, i))
 
     # Obtain corrected measurements by shaping back into a vector
     return prob_vect_t.reshape(2**N)
@@ -82,13 +107,58 @@ def renormalize_counts(corrected_counts: npt.NDArray, n_shots: int) -> npt.NDArr
 
         renormalization_factor = n_shots / sum_corrected_counts
         corrected_counts = np.rint(corrected_counts * renormalization_factor).astype(int)
+
+    # At this point, the count should be off by at most 2, added or substracted to/from the
+    # max count.
+    if sum(corrected_counts) != n_shots:
+        count_diff = sum(corrected_counts) - n_shots
+        max_count_bs = np.argmax(corrected_counts)
+        corrected_counts[max_count_bs] -= count_diff
+
     return corrected_counts
 
 
 def matrix_inv(K: npt.NDArray) -> npt.NDArray:
     return inv(K) if matrix_rank(K) == K.shape[0] else pinv(K)
 
 
+def constrained_inversion(noise_matrices: npt.NDArray, p_raw: npt.NDArray) -> npt.NDArray:
+    # Initial random guess in [0,1].
+    p_corr0 = np.random.rand(len(p_raw))
+    # Stochasticity constraints.
+    normality_constraint = LinearConstraint(np.ones(len(p_raw)).astype(int), lb=1.0, ub=1.0)
+    positivity_constraint = LinearConstraint(np.eye(len(p_raw)).astype(int), lb=0.0, ub=1.0)
+    constraints = [normality_constraint, positivity_constraint]
+    # Minimize the corrected probabilities.
+    res = minimize(corrected_probas, p_corr0, args=(noise_matrices, p_raw), constraints=constraints)
+
+    return res.x
+
+
+def majority_vote(noise_matrices: npt.NDArray, p_raw: npt.NDArray) -> npt.NDArray:
+    n_qubits = len(noise_matrices)
+    output = 0
+
+    for i in range(n_qubits):
+        temp = p_raw.reshape([2] * n_qubits)
+        transpose_axes = [i] + list(range(0, i)) + list(range(i + 1, n_qubits))
+        temp = np.transpose(temp, axes=transpose_axes).reshape(2, 2**n_qubits // 2)
+        probs = np.sum(temp, axis=1)
+
+        # given the output to be 0, the probability of observed outcomes
+        prob_zero = noise_matrices[i][1][0] ** probs[1] * noise_matrices[i][0][0] ** probs[0]
+        # given the output to be 0, the probability of observed outcomes
+        prob_one = noise_matrices[i][1][1] ** probs[1] * noise_matrices[i][0][1] ** probs[0]
+
+        if prob_one > prob_zero:
+            output += 2 ** (n_qubits - 1 - i)
+
+    p_corr = np.zeros(2**n_qubits)
+    p_corr[output] = 1
+
+    return p_corr
+
+
 def mitigation_minimization(
     noise: Noise,
     options: dict,
@@ -98,12 +168,18 @@ def mitigation_minimization(
 
     See Equation (5) in https://arxiv.org/pdf/2001.09980.pdf.
     See Page(3) in https://arxiv.org/pdf/1106.5458.pdf for MLE implementation
+    See Equation (5) in https://arxiv.org/pdf/2108.12518.pdf for MTHREE implementation
+    (matrix free measurement mitigation)
+    This method is reserved for implementations of over 20 qubits
+    See Equation (6) and Page 4 in https://arxiv.org/pdf/2402.11830.pdf for MV implementation
+    This method is reserved for algorithms with a single bit string output
+    See page (5) for extension to more than one solution (inefficient)
 
     Args:
         noise: Specifies confusion matrix and default error probability
         mitigation: Selects additional mitigation options based on noise choice.
         For readout we have the following mitigation options for optimization
-        1. constrained 2. mle. Default : mle
+        1. constrained 2. mle (Default) mle 3. mthree 4. majority vote (mv)
         samples: List of samples to be mitigated
 
     Returns:
@@ -116,48 +192,42 @@ def mitigation_minimization(
     corrected_counters: list[Counter] = []
 
     for sample in samples:
-        bitstring_length = 2**n_qubits
-        # List of bitstrings in lexicographical order.
-        ordered_bitstrings = [f"{i:0{n_qubits}b}" for i in range(bitstring_length)]
+        ordered_bitstrings = [bin(k)[2:].zfill(n_qubits) for k in range(2**n_qubits)]
         # Array of raw probabilites.
         p_raw = np.array([sample[bs] for bs in ordered_bitstrings]) / n_shots
 
         if optimization_type == ReadOutOptimization.CONSTRAINED:
-            # Initial random guess in [0,1].
-            p_corr0 = np.random.rand(bitstring_length)
-            # Stochasticity constraints.
-            normality_constraint = LinearConstraint(
-                np.ones(bitstring_length).astype(int), lb=1.0, ub=1.0
-            )
-            positivity_constraint = LinearConstraint(
-                np.eye(bitstring_length).astype(int), lb=0.0, ub=1.0
-            )
-            constraints = [normality_constraint, positivity_constraint]
-            # Minimize the corrected probabilities.
-            res = minimize(
-                corrected_probas, p_corr0, args=(noise_matrices, p_raw), constraints=constraints
-            )
-            p_corr = res.x
+            p_corr = constrained_inversion(noise_matrices, p_raw)
 
         elif optimization_type == ReadOutOptimization.MLE:
             noise_matrices_inv = list(map(matrix_inv, noise_matrices))
             # Compute corrected inverse using matrix inversion and run MLE.
             p_corr = mle_solve(tensor_rank_mult(noise_matrices_inv, p_raw))
+
+        elif optimization_type == ReadOutOptimization.MTHREE:
+            confusion_matrix_subspace = normalized_subspace_kron(noise_matrices, p_raw.nonzero()[0])
+            # GMRES (Generalized minimal residual) for linear equations in higher dimension
+            p_corr, exit_code = gmres(confusion_matrix_subspace, p_raw)
+
+            if exit_code != 0:
+                logger.warning(f"gmres did not converge, exited with code {exit_code}")
+
+            # To ensure that we are not working with negative probabilities
+            p_corr = mle_solve(p_corr)
+
+        elif optimization_type == ReadOutOptimization.MV:
+            # Majority vote : lets return just that one bit string with all the counts for now
+            p_corr = majority_vote(noise_matrices, p_raw)
+
         else:
             raise NotImplementedError(
                 f"Requested method {optimization_type} does not match supported protocols."
             )
 
         corrected_counts = np.rint(p_corr * n_shots).astype(int)
-
         # Renormalize if total counts differs from n_shots.
         corrected_counts = renormalize_counts(corrected_counts=corrected_counts, n_shots=n_shots)
-        # At this point, the count should be off by at most 2, added or substracted to/from the
-        # max count.
-        if sum(corrected_counts) != n_shots:
-            count_diff = sum(corrected_counts) - n_shots
-            max_count_bs = np.argmax(corrected_counts)
-            corrected_counts[max_count_bs] -= count_diff
+
         assert (
             corrected_counts.sum() == n_shots
         ), f"Corrected counts sum: {corrected_counts.sum()}, n_shots: {n_shots}"

qadence_protocols/types.py

@@ -25,3 +25,17 @@ class Protocols(StrEnum):
 
 
 qadence_available_protocols = Protocols.list()
+
+
+class ReadOutOptimization(StrEnum):
+    # basic inversion and maximum likelihood estimate
+    MLE = "mle"
+
+    # constrained inverse optimization
+    CONSTRAINED = "constrained"
+
+    # matrix free measurement mitigation
+    MTHREE = "mthree"
+
+    # majority voting
+    MV = "majority_vote"