[Draft] Adding Subspace Expansion error mitigation technique

mitiq/subspace_expansion/__init__.py

@@ -0,0 +1,35 @@
+# Copyright (C) 2021 Unitary Fund
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+"""Subspace Expansion method for error mitigation as introduced in:
+
+[1] McClean, J.R., Jiang, Z., Rubin, N.C., Babbush, R. and Neven, H., 2020.
+    "Decoding quantum errors with subspace expansions". 
+    (https://arxiv.org/abs/1903.05786)
+[2] Yoshioka, N., Hakoshima, H., Matsuzaki, Y., Tokunaga, Y., Suzuki, Y. and Endo, S., 2022.
+    "Generalized quantum subspace expansion".
+    (https://arxiv.org/abs/2107.02611)
+
+"""
+
+from mitiq.subspace_expansion.subspace_expansion import (
+    execute_with_subspace_expansion,
+)
+
+from mitiq.subspace_expansion.utils import (
+    convert_from_Mitiq_Observable_to_cirq_PauliSum,
+    convert_from_cirq_PauliSum_to_Mitiq_Observable,
+    convert_from_cirq_PauliString_to_Mitiq_PauliString,
+)

mitiq/subspace_expansion/subspace_expansion.py

@@ -0,0 +1,136 @@
+"""API for using Subspace Expansion error mitigation."""
+
+from typing import Callable, Sequence, Union
+from mitiq import Executor, Observable, QPROGRAM, QuantumResult, PauliString
+import cirq
+from scipy.linalg import eig
+
+from mitiq.subspace_expansion.utils import (
+    convert_from_cirq_PauliString_to_Mitiq_PauliString,
+    convert_from_Mitiq_Observable_to_cirq_PauliSum,
+)
+
+
+cache = {}  # tuple of pauli mask: expectation value
+
+
+def execute_with_subspace_expansion(
+    circuit: QPROGRAM,
+    executor: Union[Executor, Callable[[QPROGRAM], QuantumResult]],
+    check_operators: Sequence[PauliString],
+    code_hamiltonian: Observable,
+    observable: Observable,
+) -> QuantumResult:
+    """Function for the calculation of an observable from some circuit of
+        interest to be mitigated with Subspace Expansion
+    Args:
+        circuit: Quantum program to execute with error mitigation.
+        executor: Executes a circuit and returns a `QuantumResult`.
+        check_operators: List of check operators that define the stabilizer code space.
+        code_hamiltonian: Hamiltonian of the code space.
+        observable: Observable to compute the mitigated expectation value of.
+    """
+    cache.clear()
+    check_operators_cirq = [
+        check_operator._pauli for check_operator in check_operators
+    ]
+    code_hamiltonian_cirq = convert_from_Mitiq_Observable_to_cirq_PauliSum(
+        code_hamiltonian
+    )
+    observable_cirq = convert_from_Mitiq_Observable_to_cirq_PauliSum(
+        observable
+    )
+    P = get_projector(
+        circuit, executor, check_operators_cirq, code_hamiltonian_cirq
+    )
+    pop = get_expectation_value_for_observable(
+        circuit, executor, P * observable_cirq * P
+    )
+    pp = get_expectation_value_for_observable(circuit, executor, P * P)
+    return pop / pp
+
+
+def get_expectation_value_for_observable(
+    qc: QPROGRAM,
+    executor: Union[Executor, Callable[[QPROGRAM], QuantumResult]],
+    observable: Union[cirq.PauliString, cirq.PauliSum],
+):
+    def get_expectation_value_for_one_pauli(pauli_string: cirq.PauliString):
+        cache_key = tuple(pauli_string.gate.pauli_mask.tolist())
+        if cache_key == ():  # if all I's
+            return 1 * pauli_string.coefficient
+        if cache_key in cache:
+            return cache[cache_key] * pauli_string.coefficient
+        mitiq_pauli_string = (
+            convert_from_cirq_PauliString_to_Mitiq_PauliString(pauli_string)
+        )
+        return Observable(mitiq_pauli_string).expectation(qc, executor)
+
+    total_expectation_value_for_observable = 0
+
+    if type(observable) == cirq.PauliString:
+        pauli_string = observable
+        total_expectation_value_for_observable += (
+            get_expectation_value_for_one_pauli(pauli_string)
+        )
+    elif type(observable) == cirq.PauliSum:
+        pauli_sum = observable
+        for pauli_string in pauli_sum:
+            total_expectation_value_for_observable += (
+                get_expectation_value_for_one_pauli(pauli_string)
+            )
+    return total_expectation_value_for_observable
+
+
+def get_projector(qc: QPROGRAM, executor, check_operators, code_hamiltonian):
+    S = compute_overlap_matrix(qc, executor, check_operators)
+    H = compute_code_hamiltonian_overlap_matrix(qc, executor, check_operators, code_hamiltonian)
+    E, C = eig(H, S)
+    minpos = list(E).index(min(E))
+    Cs = C[:, minpos]
+
+    projector_formed_out_of_check_operators_and_Cs = sum(
+        [check_operators[i] * Cs[i] for i in range(len(check_operators))]
+    )
+
+    return projector_formed_out_of_check_operators_and_Cs
+
+
+def compute_overlap_matrix(
+    qc: QPROGRAM,
+    executor: Union[Executor, Callable[[QPROGRAM], QuantumResult]],
+    check_operators: Sequence[cirq.PauliString],
+):
+    S = []
+    # S_ij = <Ψ|Mi† Mj|Ψ>
+    for i in range(len(check_operators)):
+        S.append([])
+        for j in range(len(check_operators)):
+            paulis_and_weights = check_operators[i] * check_operators[j]
+            sij = get_expectation_value_for_observable(
+                qc, executor, paulis_and_weights
+            )
+            S[-1].append(sij)
+    return S
+
+
+def compute_code_hamiltonian_overlap_matrix(
+    qc: QPROGRAM,
+    executor: Union[Executor, Callable[[QPROGRAM], QuantumResult]],
+    check_operators: Sequence[cirq.PauliString],
+    code_hamiltonian,
+):
+    H = []
+    # H_ij = <Ψ|Mi† H Mj|Ψ>
+    for i in range(len(check_operators)):
+        H.append([])
+        for j in range(len(check_operators)):
+            paulis_and_weights = (
+                check_operators[i] * code_hamiltonian * check_operators[j]
+            )
+            H[-1].append(
+                get_expectation_value_for_observable(
+                    qc, executor, paulis_and_weights
+                )
+            )
+    return H

mitiq/subspace_expansion/tests/test_subspace_expansion.py

@@ -0,0 +1,197 @@
+# Copyright (C) 2021 Unitary Fund
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+"""Tests for the Clifford data regression top-level API."""
+import pytest
+
+from typing import List
+
+import numpy as np
+
+import cirq
+from cirq import X, Z, I
+
+from mitiq import PauliString, Observable, QPROGRAM
+
+from mitiq import SUPPORTED_PROGRAM_TYPES
+
+from mitiq.subspace_expansion import (
+    execute_with_subspace_expansion,
+)
+from mitiq.subspace_expansion.utils import (
+    convert_from_cirq_PauliSum_to_Mitiq_Observable,
+)
+
+from mitiq.interface import convert_from_mitiq, convert_to_mitiq
+
+from mitiq.interface.mitiq_cirq import compute_density_matrix
+
+
+# Allow execution with any QPROGRAM for testing.
+def execute(circuit: QPROGRAM) -> np.ndarray:
+    return compute_density_matrix(convert_to_mitiq(circuit)[0])
+
+
+def batched_execute(circuits) -> List[np.ndarray]:
+    return [execute(circuit) for circuit in circuits]
+
+
+def simulate(circuit: QPROGRAM) -> np.ndarray:
+    return compute_density_matrix(
+        convert_to_mitiq(circuit)[0], noise_level=(0,)
+    )
+
+
+@pytest.mark.parametrize("circuit_type", SUPPORTED_PROGRAM_TYPES.keys())
+def test_execute_with_subspace_expansion(circuit_type):
+    qc, actual_qubits = prepare_logical_0_state_for_5_1_3_code()
+    qc = convert_from_mitiq(qc, circuit_type)
+    (
+        check_operators,
+        code_hamiltonian,
+    ) = get_check_operators_and_code_hamiltonian()
+    observable = convert_from_cirq_PauliSum_to_Mitiq_Observable(
+        get_observable_cirq_in_code_space(actual_qubits, [Z, Z, Z, Z, Z])
+    )
+    mitigated_value = execute_with_subspace_expansion(
+        qc, execute, check_operators, code_hamiltonian, observable
+    )
+    assert abs(mitigated_value.real - 1) < 0.001
+
+    observable = convert_from_cirq_PauliSum_to_Mitiq_Observable(
+        get_observable_cirq_in_code_space(actual_qubits, [X, X, X, X, X])
+    )
+    mitigated_value = execute_with_subspace_expansion(
+        qc, execute, check_operators, code_hamiltonian, observable
+    )
+    assert abs(mitigated_value.real - 0.5) < 0.001
+
+
+def get_observable_cirq_in_code_space(
+    actual_qubits, observable: list[cirq.PauliString]
+):
+    FIVE_I = cirq.PauliString(
+        [I(actual_qubits[i]) for i in range(len(actual_qubits))]
+    )
+    projector_onto_code_space = [
+        [X, Z, Z, X, I],
+        [I, X, Z, Z, X],
+        [X, I, X, Z, Z],
+        [Z, X, I, X, Z],
+    ]
+
+    observable_in_code_space = FIVE_I
+    all_paulis = projector_onto_code_space + [observable]
+    for g in all_paulis:
+        f = cirq.PauliString([(g[i])(actual_qubits[i]) for i in range(len(g))])
+        observable_in_code_space *= 0.5 * (FIVE_I + f)
+    return observable_in_code_space
+
+
+def get_check_operators_and_code_hamiltonian() -> tuple:
+    Ms = [
+        "YIYXX",
+        "ZIZYY",
+        "IXZZX",
+        "ZXIXZ",
+        "YYZIZ",
+        "XYIYX",
+        "YZIZY",
+        "ZZXIX",
+        "XZZXI",
+        "ZYYZI",
+        "IYXXY",
+        "IZYYZ",
+        "YXXYI",
+        "XXYIY",
+        "XIXZZ",
+        "IIIII",
+    ]
+    Ms_as_pauliStrings = [
+        PauliString(M, coeff=1, support=range(5)) for M in Ms
+    ]
+    negative_Ms_as_pauliStrings = [
+        PauliString(M, coeff=-1, support=range(5)) for M in Ms
+    ]
+    Hc = Observable(*negative_Ms_as_pauliStrings)
+    return Ms_as_pauliStrings, Hc
+
+
+def prepare_logical_0_state_for_5_1_3_code():
+    def gram_schmidt(
+        orthogonal_vecs: List[np.ndarray],
+    ):
+        # normalize input
+        orthonormalVecs = [
+            vec / np.sqrt(np.vdot(vec, vec)) for vec in orthogonal_vecs
+        ]
+        dim = np.shape(orthogonal_vecs[0])[0]  # get dim of vector space
+        for i in range(dim - len(orthogonal_vecs)):
+            new_vec = np.zeros(dim)
+            new_vec[i] = 1  # construct ith basis vector
+            projs = sum(
+                [
+                    np.vdot(new_vec, cached_vec) * cached_vec
+                    for cached_vec in orthonormalVecs
+                ]
+            )  # sum of projections of new vec with all existing vecs
+            new_vec -= projs
+            orthonormalVecs.append(
+                new_vec / np.sqrt(np.vdot(new_vec, new_vec))
+            )
+        return orthonormalVecs
+
+    logical_0_state = np.zeros(32)
+    for z in ["00000", "10010", "01001", "10100", "01010", "00101"]:
+        logical_0_state[int(z, 2)] = 1 / 4
+    for z in [
+        "11011",
+        "00110",
+        "11000",
+        "11101",
+        "00011",
+        "11110",
+        "01111",
+        "10001",
+        "01100",
+        "10111",
+    ]:
+        logical_0_state[int(z, 2)] = -1 / 4
+
+    logical_1_state = np.zeros(32)
+    for z in ["11111", "01101", "10110", "01011", "10101", "11010"]:
+        logical_1_state[int(z, 2)] = 1 / 4
+    for z in [
+        "00100",
+        "11001",
+        "00111",
+        "00010",
+        "11100",
+        "00001",
+        "10000",
+        "01110",
+        "10011",
+        "01000",
+    ]:
+        logical_1_state[int(z, 2)] = -1 / 4
+
+    # Fill up the rest of the matrix with orthonormal vectors
+    res = gram_schmidt([logical_0_state, logical_1_state])
+    mat = np.reshape(res, (32, 32)).T
+    circuit = cirq.Circuit()
+    g = cirq.MatrixGate(mat)
+    actual_qubits = cirq.LineQubit.range(5)
+    circuit.append(g(*actual_qubits))
+    return circuit, actual_qubits