diff --git a/doc/source/api-reference/qibo.rst b/doc/source/api-reference/qibo.rst index 669fe34c6c..bbeeda603c 100644 --- a/doc/source/api-reference/qibo.rst +++ b/doc/source/api-reference/qibo.rst @@ -279,7 +279,7 @@ and both can be used as standalone functions or in combination with the other general mitigation methods by setting the paramter `readout`. -Calibration Matrix +Response Matrix """""""""""""""""" Given :math:`n` qubits, all the possible :math:`2^n` states are constructed via the application of the corresponding sequence of :math:`X` gates @@ -288,7 +288,7 @@ In the presence of readout errors, we will measure for each state :math:`i` some frequencies :math:`F_i^{noisy}` different from the ideal ones :math:`F_i^{ideal}=\delta_{i,j}`. -The effect of the error is modeled by the matrix composed of the noisy frequencies as +The effect of the error is modeled by the response matrix composed of the noisy frequencies as columns :math:`M=\big(F_0^{noisy},...,F_{n-1}^{noisy}\big)`. We have indeed that: .. math:: @@ -297,14 +297,31 @@ columns :math:`M=\big(F_0^{noisy},...,F_{n-1}^{noisy}\big)`. We have indeed that and, therefore, the calibration matrix obtained as :math:`M_{\text{cal}}=M^{-1}` can be used to recover the noise-free frequencies. -.. autofunction:: qibo.models.error_mitigation.calibration_matrix +The calibration matrix :math:`M_{\text{cal}}` lacks stochasticity, resulting in a 'negative probability' issue. +The distributions that arise after applying :math:`M_{\text{cal}}` are quasiprobabilities; +the individual elements can be negative surpass 1, provided they sum to 1. +It is posible to use Iterative Bayesian Unfolding (IBU) to preserve non-negativity. +See `Nachman et al `_ for more details. -.. autofunction:: qibo.models.error_mitigation.apply_readout_mitigation +.. autofunction:: qibo.models.error_mitigation.get_response_matrix -Randomized -"""""""""" + +.. autofunction:: qibo.models.error_mitigation.iterative_bayesian_unfolding + + +.. autofunction:: qibo.models.error_mitigation.apply_resp_mat_readout_mitigation + + +.. autofunction:: qibo.models.error_mitigation.apply_randomized_readout_mitigation + + +.. autofunction:: qibo.models.error_mitigation.get_expectation_val_with_readout_mitigation + + +Randomized readout mitigation +"""""""""""""""""""""""""""""" This approach converts the effect of any noise map :math:`A` into a single multiplication factor for each Pauli observable, that is, diagonalizes the measurement channel. The multiplication factor :math:`\lambda` can be directly measured even without @@ -320,7 +337,7 @@ factor results in the mitigated Pauli expectation value :math:`\langle O\rangle_ Zero Noise Extrapolation (ZNE) """""""""""""""""""""""""""""" -Given a noisy circuit :math:`C` and an observable :math:`A`, Zero Noise Extrapolation (ZNE) +Given a noisy circuit :math:`C` and an observable :math:`A`, Zero Noise Extrapolation (ZNE) consists in running :math:`n+1` versions of the circuit with different noise levels :math:`\{c_j\}_{j=0..n}` and, for each of them, measuring the expected value of the observable :math:`E_j=\langle A\rangle_j`. @@ -381,7 +398,7 @@ See `Sopena et al `_ for more details. .. autofunction:: qibo.models.error_mitigation.CDR -.. autofunction:: qibo.models.error_mitigation.sample_training_circuit +.. autofunction:: qibo.models.error_mitigation.sample_training_circuit_cdr Variable Noise CDR (vnCDR) @@ -416,6 +433,34 @@ See `Sopena et al `_ for all the details. .. autofunction:: qibo.models.error_mitigation.vnCDR +Importance Clifford Sampling (ICS) +"""""""""""""""""""""""""""""""""" + +In the Importance Clifford Sampling (ICS) method, a set of :math:`n` circuits +:math:`S_n=\{C_i\}_{i=1,..,n}` that stabilizes a given Pauli observable is generated starting from the original circuit +:math:`C_0` by replacing all the non-Clifford gates with Clifford ones. +Given an observable :math:`A`, all the circuits of :math:`S_n` are both simulated +to obtain the correspondent expected values of :math:`A` in noise-free condition +:math:`\{a_i^{exact}\}_{i=1,..,n}`, and run in noisy conditions to obtain the noisy +expected values :math:`\{a_i^{noisy}\}_{i=1,..,n}`. + +Finally, a theoretically inspired model :math:`f` is learned using the training data. + +The mitigated expected value of :math:`A` at the end of :math:`C_0` is then +obtained simply with :math:`f(a_0^{noisy})`. + +In this implementation the initial circuit is expected to be decomposed in the three +Clifford gates :math:`RX(\frac{\pi}{2})`, :math:`CNOT`, :math:`X` and in :math:`RZ(\theta)` +(which is Clifford only for :math:`\theta=\frac{n\pi}{2}`). +By default the set of Clifford gates used for substitution is +:math:`\{RZ(0),RZ(\frac{\pi}{2}),RZ(\pi),RZ(\frac{3}{2}\pi)\}`. +See `Sopena et al `_ for more details. + +.. autofunction:: qibo.models.error_mitigation.ICS + + +.. autofunction:: qibo.models.error_mitigation.sample_clifford_training_circuit + _______________________ .. _Gates: @@ -1420,10 +1465,6 @@ passing a symplectic matrix to the constructor. symplectic_matrix = backend.zero_state(nqubits=3) clifford = Clifford(symplectic_matrix, engine=NumpyBackend()) - # The initialization above is equivalent to the initialization below - circuit = Circuit(nqubits=3) - clifford = Clifford(circuit, engine=NumpyBackend()) - The generators of the stabilizers can be extracted with the :meth:`qibo.quantum_info.clifford.Clifford.generators` method, or the complete set of :math:`d = 2^{n}` stabilizers operators can be extracted through the diff --git a/doc/source/code-examples/advancedexamples.rst b/doc/source/code-examples/advancedexamples.rst index de00889d14..ca85dd2e25 100644 --- a/doc/source/code-examples/advancedexamples.rst +++ b/doc/source/code-examples/advancedexamples.rst @@ -1203,23 +1203,23 @@ Let's see how to use them. For starters, let's define a dummy circuit with some hz = 0.5 hx = 0.5 dt = 0.25 - c = Circuit(nqubits, density_matrix=True) - c.add(gates.RZ(q, theta=-2 * hz * dt - np.pi / 2) for q in range(nqubits)) - c.add(gates.RX(q, theta=np.pi / 2) for q in range(nqubits)) - c.add(gates.RZ(q, theta=-2 * hx * dt + np.pi) for q in range(nqubits)) - c.add(gates.RX(q, theta=np.pi / 2) for q in range(nqubits)) - c.add(gates.RZ(q, theta=-np.pi / 2) for q in range(nqubits)) - c.add(gates.CNOT(q, q + 1) for q in range(0, nqubits - 1, 2)) - c.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(0, nqubits - 1, 2)) - c.add(gates.CNOT(q, q + 1) for q in range(0, nqubits - 1, 2)) - c.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2)) - c.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(1, nqubits, 2)) - c.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2)) + circ = Circuit(nqubits, density_matrix=True) + circ.add(gates.RZ(q, theta=-2 * hz * dt - np.pi / 2) for q in range(nqubits)) + circ.add(gates.RX(q, theta=np.pi / 2) for q in range(nqubits)) + circ.add(gates.RZ(q, theta=-2 * hx * dt + np.pi) for q in range(nqubits)) + circ.add(gates.RX(q, theta=np.pi / 2) for q in range(nqubits)) + circ.add(gates.RZ(q, theta=-np.pi / 2) for q in range(nqubits)) + circ.add(gates.CNOT(q, q + 1) for q in range(0, nqubits - 1, 2)) + circ.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(0, nqubits - 1, 2)) + circ.add(gates.CNOT(q, q + 1) for q in range(0, nqubits - 1, 2)) + circ.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2)) + circ.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(1, nqubits, 2)) + circ.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2)) # Include the measurements - c.add(gates.M(*range(nqubits))) + circ.add(gates.M(*range(nqubits))) # visualize the circuit - print(c.draw()) + print(circ.draw()) # q0: ─RZ─RX─RZ─RX─RZ─o────o────────M─ # q1: ─RZ─RX─RZ─RX─RZ─X─RZ─X─o────o─M─ @@ -1252,7 +1252,7 @@ the real quantum hardware, instead, we can use a noise model: .. testcode:: # Noise-free expected value - exact = obs.expectation(backend.execute_circuit(c).state()) + exact = obs.expectation(backend.execute_circuit(circ).state()) print(exact) # 0.9096065335014379 @@ -1270,7 +1270,7 @@ the real quantum hardware, instead, we can use a noise model: ) noise.add(ReadoutError(probabilities=prob), gate=gates.M) # Noisy expected value without mitigation - noisy = obs.expectation(backend.execute_circuit(noise.apply(c)).state()) + noisy = obs.expectation(backend.execute_circuit(noise.apply(circ)).state()) print(noisy) # 0.5647937721701448 @@ -1287,22 +1287,21 @@ Now let's check that error mitigation produces better estimates of the exact exp Readout Mitigation ^^^^^^^^^^^^^^^^^^ Firstly, let's try to mitigate the readout errors. To do this, we can either compute the -calibration matrix and use it modify the final state after the circuit execution: +response matrix and use it modify the final state after the circuit execution: .. testcode:: - from qibo.models.error_mitigation import apply_readout_mitigation, calibration_matrix + from qibo.models.error_mitigation import get_expectation_val_with_readout_mitigation, get_response_matrix nshots = 10000 - # compute the calibration matrix - calibration = calibration_matrix( + # compute the response matrix + response_matrix = get_response_matrix( nqubits, backend=backend, noise_model=noise, nshots=nshots ) - # execute the circuit - state = backend.execute_circuit(noise.apply(c), nshots=nshots) + # define mitigation options + readout = {"response_matrix": response_matrix} # mitigate the readout errors - mit_state = apply_readout_mitigation(state, calibration) - mit_val = mit_state.expectation_from_samples(obs) + mit_val = get_expectation_val_with_readout_mitigation(circ, obs, noise, readout=readout) print(mit_val) # 0.5945794816381054 @@ -1317,13 +1316,10 @@ Or use the randomized readout mitigation: from qibo.models.error_mitigation import apply_randomized_readout_mitigation - ncircuits = 10 - result, result_cal = apply_randomized_readout_mitigation( - c, backend=backend, noise_model=noise, nshots=nshots, ncircuits=ncircuits - ) - mit_val = result.expectation_from_samples( - obs - ) / result_cal.expectation_from_samples(obs) + # define mitigation options + readout = {"ncircuits": 10} + # mitigate the readout errors + mit_val = get_expectation_val_with_readout_mitigation(circ, obs, noise, readout=readout) print(mit_val) # 0.5860884499785314 @@ -1363,7 +1359,7 @@ For example if we use the five levels ``[0,1,2,3,4]`` : # Mitigated expected value estimate = ZNE( - circuit=c, + circuit=circ, observable=obs, noise_levels=np.arange(5), noise_model=noise, @@ -1385,14 +1381,14 @@ combined with the readout mitigation: .. testcode:: # we can either use - # the calibration matrix computed earlier - readout = {'calibration_matrix': calibration} + # the response matrix computed earlier + readout = {'response_matrix': response_matrix} # or the randomized readout readout = {'ncircuits': 10} # Mitigated expected value estimate = ZNE( - circuit=c, + circuit=circ, observable=obs, backend=backend, noise_levels=np.arange(5), @@ -1422,15 +1418,16 @@ circuit is expected to be decomposed in the set of primitive gates :math:`RX(\fr # Mitigated expected value estimate = CDR( - circuit=c, + circuit=circ, observable=obs, + n_training_samples=10, backend=backend, noise_model=noise, nshots=10000, readout=readout, ) print(estimate) - # 0.9090604794014961 + # 0.8983676333969615 .. testoutput:: :hide: @@ -1439,6 +1436,7 @@ circuit is expected to be decomposed in the set of primitive gates :math:`RX(\fr Again, the mitigated expected value improves over the noisy one and is also slightly better compared to ZNE. + Variable Noise CDR (vnCDR) ^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -1451,8 +1449,9 @@ caveat about the input circuit for CDR is valid here as well. # Mitigated expected value estimate = vnCDR( - circuit=c, + circuit=circ, observable=obs, + n_training_samples=10, backend=backend, noise_levels=np.arange(3), noise_model=noise, @@ -1461,7 +1460,7 @@ caveat about the input circuit for CDR is valid here as well. readout=readout, ) print(estimate) - # 0.9085991439303123 + # 0.8998376314644383 .. testoutput:: :hide: @@ -1472,6 +1471,35 @@ The result is similar to the one obtained by CDR. Usually, one would expect slig however, this can substantially vary depending on the circuit and the observable considered and, therefore, it is hard to tell a priori. + +Importance Clifford Sampling (ICS) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The use of iCS is straightforward, analogous to CDR and vnCDR. + +.. testcode:: + + from qibo.models.error_mitigation import ICS + + # Mitigated expected value + estimate = ICS( + circuit=circ, + observable=obs, + n_training_samples=10, + backend=backend, + noise_model=noise, + nshots=10000, + readout=readout, + ) + print(estimate) + # 0.9183495097398502 + +.. testoutput:: + :hide: + + ... + +Again, the mitigated expected value improves over the noisy one and is also slightly better compared to ZNE. This was just a basic example usage of the three methods, for all the details about them you should check the API-reference page :ref:`Error Mitigation `. .. _timeevol-example: diff --git a/src/qibo/models/__init__.py b/src/qibo/models/__init__.py index ad0b48bbde..cf55bf7d52 100644 --- a/src/qibo/models/__init__.py +++ b/src/qibo/models/__init__.py @@ -1,14 +1,7 @@ from qibo.models import hep, tsp from qibo.models.circuit import Circuit from qibo.models.encodings import unary_encoder -from qibo.models.error_mitigation import ( - CDR, - ZNE, - get_gammas, - get_noisy_circuit, - sample_training_circuit, - vnCDR, -) +from qibo.models.error_mitigation import CDR, ICS, ZNE, vnCDR from qibo.models.evolution import AdiabaticEvolution, StateEvolution from qibo.models.grover import Grover from qibo.models.qft import QFT diff --git a/src/qibo/models/error_mitigation.py b/src/qibo/models/error_mitigation.py index e29619c3e0..2e2d577d2e 100644 --- a/src/qibo/models/error_mitigation.py +++ b/src/qibo/models/error_mitigation.py @@ -1,6 +1,6 @@ """Error Mitigation Methods.""" -from math import factorial +from itertools import product import numpy as np from scipy.optimize import curve_fit @@ -10,42 +10,45 @@ from qibo.config import raise_error -def get_gammas(c, solve: bool = True): +def get_gammas(noise_levels, analytical: bool = True): """Standalone function to compute the ZNE coefficients given the noise levels. Args: - c (numpy.ndarray): array containing the different noise levels. + noise_levels (numpy.ndarray): array containing the different noise levels. Note that in the CNOT insertion paradigm this corresponds to the number of CNOT pairs to be inserted. The canonical ZNE noise levels are obtained as ``2 * c + 1``. - solve (bool, optional): If ``True``, computes the coeffients by solving the + analytical (bool, optional): if ``True``, computes the coeffients by solving the linear system. If ``False``, use the analytical solution valid for the CNOT insertion method. Default is ``True``. Returns: - numpy.ndarray: The computed coefficients. + numpy.ndarray: the computed coefficients. """ - if solve: - c = 2 * c + 1 - a = np.array([c**i for i in range(len(c))]) - b = np.zeros(len(c)) - b[0] = 1 - gammas = np.linalg.solve(a, b) + if analytical: + noise_levels = 2 * noise_levels + 1 + a_matrix = np.array([noise_levels**i for i in range(len(noise_levels))]) + b_vector = np.zeros(len(noise_levels)) + b_vector[0] = 1 + zne_coefficients = np.linalg.solve(a_matrix, b_vector) else: - cmax = c[-1] - gammas = np.array( + max_noise_level = noise_levels[-1] + zne_coefficients = np.array( [ 1 - / (2 ** (2 * cmax) * factorial(i)) + / (2 ** (2 * max_noise_level) * np.math.factorial(i)) * (-1) ** i / (1 + 2 * i) - * factorial(1 + 2 * cmax) - / (factorial(cmax) * factorial(cmax - i)) - for i in c + * np.math.factorial(1 + 2 * max_noise_level) + / ( + np.math.factorial(max_noise_level) + * np.math.factorial(max_noise_level - i) + ) + for i in noise_levels ] ) - return gammas + return zne_coefficients def get_noisy_circuit(circuit, num_insertions: int, insertion_gate: str = "CNOT"): @@ -59,7 +62,7 @@ def get_noisy_circuit(circuit, num_insertions: int, insertion_gate: str = "CNOT" Default is ``"CNOT"``. Returns: - :class:`qibo.models.Circuit`: The circuit with the inserted CNOT pairs. + :class:`qibo.models.Circuit`: circuit with the inserted gate pairs. """ if insertion_gate not in ("CNOT", "RX"): # pragma: no cover raise_error( @@ -80,16 +83,17 @@ def get_noisy_circuit(circuit, num_insertions: int, insertion_gate: str = "CNOT" for gate in circuit.queue: noisy_circuit.add(gate) + if isinstance(gate, i_gate): if insertion_gate == "CNOT": control = gate.control_qubits[0] target = gate.target_qubits[0] - for i in range(num_insertions): + for _ in range(num_insertions): noisy_circuit.add(gates.CNOT(control, target)) noisy_circuit.add(gates.CNOT(control, target)) elif gate.init_kwargs["theta"] == theta: qubit = gate.qubits[0] - for i in range(num_insertions): + for _ in range(num_insertions): noisy_circuit.add(gates.RX(qubit, theta=theta)) noisy_circuit.add(gates.RX(qubit, theta=-theta)) @@ -101,10 +105,11 @@ def ZNE( observable, noise_levels, noise_model=None, - nshots=int(1e4), + nshots=10000, solve_for_gammas=False, insertion_gate="CNOT", - readout: dict = {}, + readout=None, + qubit_map=None, backend=None, ): """Runs the Zero Noise Extrapolation method for error mitigation. @@ -115,57 +120,61 @@ def ZNE( Args: circuit (:class:`qibo.models.Circuit`): input circuit. - observable (numpy.ndarray): Observable to measure. + observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): Observable to measure. noise_levels (numpy.ndarray): Sequence of noise levels. noise_model (:class:`qibo.noise.NoiseModel`, optional): Noise model applied to simulate noisy computation. - nshots (int, optional): Number of shots. + nshots (int, optional): Number of shots. Defaults to :math:`10000`. solve_for_gammas (bool, optional): If ``True``, explicitly solve the - equations to obtain the ``gamma`` coefficients. + equations to obtain the ``gamma`` coefficients. Default is ``False``. insertion_gate (str, optional): gate to be used in the insertion. If ``"RX"``, the gate used is :math:``RX(\\pi / 2)``. - Default is ``"CNOT"``. - readout (dict, optional): It has the structure - {'calibration_matrix': `numpy.ndarray`, 'ncircuits': `int`}. - If passed, the calibration matrix or the randomized method is - used to mitigate readout errors. - backend (:class:`qibo.backends.abstract.Backend`, optional): Calculation engine. + Defaults to ``"CNOT"``. + readout (dict, optional): a dictionary that may contain the following keys: + + * ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation. + * response_matrix: numpy.ndarray, used for applying a pre-computed response matrix for readout error mitigation. + * ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}. + + qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``. + backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used + in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. Returns: numpy.ndarray: Estimate of the expected value of ``observable`` in the noise free condition. - """ + Reference: + 1. K. Temme, S. Bravyi et al, *Error mitigation for short-depth quantum circuits*. + `arXiv:1612.02058 [quant-ph] `_. + """ if backend is None: # pragma: no cover backend = GlobalBackend() + if readout is None: + readout = {} - expected_val = [] + expected_values = [] for num_insertions in noise_levels: noisy_circuit = get_noisy_circuit( circuit, num_insertions, insertion_gate=insertion_gate ) - if "ncircuits" in readout.keys(): - circuit_result, circuit_result_cal = apply_randomized_readout_mitigation( - noisy_circuit, noise_model, nshots, readout["ncircuits"], backend - ) - else: - if noise_model is not None and backend.name != "qibolab": - noisy_circuit = noise_model.apply(noisy_circuit) - circuit_result = backend.execute_circuit(noisy_circuit, nshots=nshots) - if "calibration_matrix" in readout.keys() is not None: - circuit_result = apply_readout_mitigation( - circuit_result, readout["calibration_matrix"] - ) - val = circuit_result.expectation_from_samples(observable) - if "ncircuits" in readout.keys(): - val /= circuit_result_cal.expectation_from_samples(observable) - expected_val.append(val) + val = get_expectation_val_with_readout_mitigation( + noisy_circuit, + observable, + noise_model, + nshots, + readout, + qubit_map, + backend=backend, + ) + expected_values.append(val) - gamma = get_gammas(noise_levels, solve=solve_for_gammas) + gamma = get_gammas(noise_levels, analytical=solve_for_gammas) - return np.sum(gamma * expected_val) + return np.sum(gamma * expected_values) -def sample_training_circuit( +def sample_training_circuit_cdr( circuit, replacement_gates: list = None, sigma: float = 0.5, @@ -181,7 +190,9 @@ def sample_training_circuit( form (``gates.XYZ``, ``kwargs``). For example, phase gates are used by default: ``list((RZ, {'theta':0}), (RZ, {'theta':pi/2}), (RZ, {'theta':pi}), (RZ, {'theta':3*pi/2}))``. sigma (float, optional): standard devation of the Gaussian distribution used for sampling. - backend (:class:`qibo.backends.abstract.Backend`, optional): Calculation engine. + backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used + in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. Returns: :class:`qibo.models.Circuit`: The sampled circuit. @@ -192,18 +203,15 @@ def sample_training_circuit( if replacement_gates is None: replacement_gates = [(gates.RZ, {"theta": n * np.pi / 2}) for n in range(4)] - # Find all the non-Clifford RZ gates gates_to_replace = [] for i, gate in enumerate(circuit.queue): if isinstance(gate, gates.RZ): if gate.init_kwargs["theta"] % (np.pi / 2) != 0.0: gates_to_replace.append((i, gate)) - if len(gates_to_replace) == 0: + if not gates_to_replace: raise_error(ValueError, "No non-Clifford RZ gate found, no circuit sampled.") - # For each RZ gate build the possible candidates and - # compute the frobenius distance to the candidates replacement, distance = [], [] for _, gate in gates_to_replace: rep_gates = np.array( @@ -221,31 +229,28 @@ def sample_training_circuit( ) distance = np.vstack(distance) - # Compute the scores prob = np.exp(-(distance**2) / sigma**2) - # Sample which of the RZ found to substitute + index = np.random.choice( range(len(gates_to_replace)), size=min(int(len(gates_to_replace) / 2), 50), replace=False, p=prob.sum(-1) / prob.sum(), ) + gates_to_replace = np.array([gates_to_replace[i] for i in index]) prob = [prob[i] for i in index] - # Sample which replacement gate to substitute with + replacement = np.array([replacement[i] for i in index]) replacement = [ replacement[i][np.random.choice(range(len(p)), size=1, p=p / p.sum())[0]] for i, p in enumerate(prob) ] replacement = {i[0]: g for i, g in zip(gates_to_replace, replacement)} - # Build the training circuit by substituting the sampled gates + sampled_circuit = circuit.__class__(**circuit.init_kwargs) for i, gate in enumerate(circuit.queue): - if i in replacement.keys(): - sampled_circuit.add(replacement[i]) - else: - sampled_circuit.add(gate) + sampled_circuit.add(replacement.get(i, gate)) return sampled_circuit @@ -254,11 +259,12 @@ def CDR( circuit, observable, noise_model, - nshots: int = int(1e4), + nshots: int = 10000, model=lambda x, a, b: a * x + b, n_training_samples: int = 100, full_output: bool = False, - readout: dict = {}, + readout=None, + qubit_map=None, backend=None, ): """Runs the Clifford Data Regression error mitigation method. @@ -266,82 +272,66 @@ def CDR( Args: circuit (:class:`qibo.models.Circuit`): input circuit decomposed in the primitive gates ``X``, ``CNOT``, ``RX(pi/2)``, ``RZ(theta)``. - observable (numpy.ndarray): observable to be measured. + observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): observable to be measured. noise_model (:class:`qibo.noise.NoiseModel`): noise model used for simulating noisy computation. - nshots (int, optional): number of shots. + nshots (int, optional): number of shots. Defaults :math:`10000`. model (callable, optional): model used for fitting. This should be a callable function object ``f(x, *params)``, taking as input the predictor variable and the parameters. Default is a simple linear model ``f(x,a,b) := a*x + b``. - n_training_samples (int, optional): number of training circuits to sample. + n_training_samples (int, optional): number of training circuits to sample. Defaults to 100. full_output (bool, optional): if ``True``, this function returns additional - information: ``val``, ``optimal_params``, ``train_val``. - readout (dict, optional): It has the structure - {'calibration_matrix': `numpy.ndarray`, 'ncircuits': `int`}. - If passed, the calibration matrix or the randomized method is - used to mitigate readout errors. - backend (:class:`qibo.backends.abstract.Backend`, optional): calculation engine. + information: ``val``, ``optimal_params``, ``train_val``. Defaults to ``False``. + readout (dict, optional): a dictionary that may contain the following keys: + + * ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation. + * response_matrix: numpy.ndarray, used for applying a pre-computed response matrix for readout error mitigation. + * ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}. + + qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``. + backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used + in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. Returns: mit_val (float): Mitigated expectation value of `observable`. val (float): Noisy expectation value of `observable`. optimal_params (list): Optimal values for `params`. train_val (dict): Contains the noise-free and noisy expectation values obtained with the training circuits. - """ - # Set backend + Reference: + 1. P. Czarnik, A. Arrasmith et al, *Error mitigation with Clifford quantum-circuit data*. + `arXiv:2005.10189 [quant-ph] `_. + """ if backend is None: # pragma: no cover backend = GlobalBackend() - # Sample the training set + if readout is None: + readout = {} + training_circuits = [ - sample_training_circuit(circuit) for n in range(n_training_samples) + sample_training_circuit_cdr(circuit) for _ in range(n_training_samples) ] - # Run the sampled circuits + train_val = {"noise-free": [], "noisy": []} - for c in training_circuits: - val = c(nshots=nshots).expectation_from_samples(observable) + for circ in training_circuits: + val = circ(nshots=nshots).expectation_from_samples(observable) train_val["noise-free"].append(val) - if "ncircuits" in readout.keys(): - circuit_result, circuit_result_cal = apply_randomized_readout_mitigation( - c, noise_model, nshots, readout["ncircuits"], backend - ) - else: - if noise_model is not None and backend.name != "qibolab": - c = noise_model.apply(c) - circuit_result = backend.execute_circuit(c, nshots=nshots) - if "calibration_matrix" in readout.keys() is not None: - circuit_result = apply_readout_mitigation( - circuit_result, readout["calibration_matrix"] - ) - val = circuit_result.expectation_from_samples(observable) - if "ncircuits" in readout.keys(): - val /= circuit_result_cal.expectation_from_samples(observable) + val = get_expectation_val_with_readout_mitigation( + circ, observable, noise_model, nshots, readout, qubit_map, backend=backend + ) train_val["noisy"].append(val) - # Fit the model + optimal_params = curve_fit(model, train_val["noisy"], train_val["noise-free"])[0] - # Run the input circuit - if "ncircuits" in readout.keys(): - circuit_result, circuit_result_cal = apply_randomized_readout_mitigation( - circuit, noise_model, nshots, readout["ncircuits"], backend - ) - else: - if noise_model is not None and backend.name != "qibolab": - circuit = noise_model.apply(circuit) - circuit_result = backend.execute_circuit(circuit, nshots=nshots) - if "calibration_matrix" in readout.keys() is not None: - circuit_result = apply_readout_mitigation( - circuit_result, readout["calibration_matrix"] - ) - val = circuit_result.expectation_from_samples(observable) - if "ncircuits" in readout.keys(): - val /= circuit_result_cal.expectation_from_samples(observable) + + val = get_expectation_val_with_readout_mitigation( + circuit, observable, noise_model, nshots, readout, qubit_map, backend=backend + ) mit_val = model(val, *optimal_params) - # Return data - if full_output == True: + if full_output is True: return mit_val, val, optimal_params, train_val - else: - return mit_val + + return mit_val def vnCDR( @@ -349,12 +339,13 @@ def vnCDR( observable, noise_levels, noise_model, - nshots: int = int(1e4), + nshots: int = 10000, model=lambda x, *params: (x * np.array(params).reshape(-1, 1)).sum(0), n_training_samples: int = 100, insertion_gate: str = "CNOT", full_output: bool = False, - readout: dict = {}, + readout=None, + qubit_map=None, backend=None, ): """Runs the variable-noise Clifford Data Regression error mitigation method. @@ -362,11 +353,11 @@ def vnCDR( Args: circuit (:class:`qibo.models.Circuit`): input circuit decomposed in the primitive gates ``X``, ``CNOT``, ``RX(pi/2)``, ``RZ(theta)``. - observable (numpy.ndarray): observable to be measured. + observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): observable to be measured. noise_levels (numpy.ndarray): sequence of noise levels. noise_model (:class:`qibo.noise.NoiseModel`): noise model used for simulating noisy computation. - nshots (int, optional): number of shots. + nshots (int, optional): number of shots. Defaults to :math:`10000`. model (callable, optional): model used for fitting. This should be a callable function object ``f(x, *params)``, taking as input the predictor variable and the parameters. Default is a simple linear model ``f(x,a,b) := a*x + b``. @@ -375,12 +366,17 @@ def vnCDR( If ``"RX"``, the gate used is :math:``RX(\\pi / 2)``. Default is ``"CNOT"``. full_output (bool, optional): if ``True``, this function returns additional - information: ``val``, ``optimal_params``, ``train_val``. - readout (dict, optional): It has the structure - {'calibration_matrix': `numpy.ndarray`, 'ncircuits': `int`}. - If passed, the calibration matrix or the randomized method is - used to mitigate readout errors. - backend (:class:`qibo.backends.abstract.Backend`, optional): calculation engine. + information: ``val``, ``optimal_params``, ``train_val``. Defaults to ``False``. + readout (dict, optional): a dictionary that may contain the following keys: + + * ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation. + * response_matrix: numpy.ndarray, used for applying a pre-computed response matrix for readout error mitigation. + * ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}. + + qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``. + backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used + in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. Returns: mit_val (float): Mitigated expectation value of `observable`. @@ -388,169 +384,213 @@ def vnCDR( optimal_params (list): Optimal values for `params`. train_val (dict): Contains the noise-free and noisy expectation values obtained with the training circuits. - """ - # Set backend + Reference: + 1. A. Lowe, MH. Gordon et al, *Unified approach to data-driven quantum error mitigation*. + `arXiv:2011.01157 [quant-ph] `_. + """ if backend is None: # pragma: no cover backend = GlobalBackend() + if readout is None: + readout = {} - # Sample the training circuits training_circuits = [ - sample_training_circuit(circuit) for n in range(n_training_samples) + sample_training_circuit_cdr(circuit) for _ in range(n_training_samples) ] train_val = {"noise-free": [], "noisy": []} - # Add the different noise levels and run the circuits - for c in training_circuits: - val = c(nshots=nshots).expectation_from_samples(observable) + for circ in training_circuits: + val = circ(nshots=nshots).expectation_from_samples(observable) train_val["noise-free"].append(val) for level in noise_levels: - noisy_c = get_noisy_circuit(c, level, insertion_gate=insertion_gate) - if "ncircuits" in readout.keys(): - ( - circuit_result, - circuit_result_cal, - ) = apply_randomized_readout_mitigation( - noisy_c, noise_model, nshots, readout["ncircuits"], backend - ) - else: - if noise_model is not None and backend.name != "qibolab": - noisy_c = noise_model.apply(noisy_c) - circuit_result = backend.execute_circuit(noisy_c, nshots=nshots) - if "calibration_matrix" in readout.keys(): - circuit_result = apply_readout_mitigation( - circuit_result, readout["calibration_matrix"] - ) - val = circuit_result.expectation_from_samples(observable) - if "ncircuits" in readout.keys(): - val /= circuit_result_cal.expectation_from_samples(observable) + noisy_c = get_noisy_circuit(circ, level, insertion_gate=insertion_gate) + val = get_expectation_val_with_readout_mitigation( + noisy_c, + observable, + noise_model, + nshots, + readout, + qubit_map, + backend=backend, + ) train_val["noisy"].append(val) - # Repeat noise-free values for each noise level noisy_array = np.array(train_val["noisy"]).reshape(-1, len(noise_levels)) - # Fit the model params = np.random.rand(len(noise_levels)) optimal_params = curve_fit(model, noisy_array.T, train_val["noise-free"], p0=params) - # Run the input circuit val = [] for level in noise_levels: noisy_c = get_noisy_circuit(circuit, level, insertion_gate=insertion_gate) - if "ncircuits" in readout.keys(): - circuit_result, circuit_result_cal = apply_randomized_readout_mitigation( - noisy_c, noise_model, nshots, readout["ncircuits"], backend - ) - else: - if noise_model is not None and backend.name != "qibolab": - noisy_c = noise_model.apply(noisy_c) - circuit_result = backend.execute_circuit(noisy_c, nshots=nshots) - if "calibration_matrix" in readout.keys(): - circuit_result = apply_readout_mitigation( - circuit_result, readout["calibration_matrix"] - ) - expval = circuit_result.expectation_from_samples(observable) - if "ncircuits" in readout.keys(): - expval /= circuit_result_cal.expectation_from_samples(observable) + expval = get_expectation_val_with_readout_mitigation( + noisy_c, + observable, + noise_model, + nshots, + readout, + qubit_map, + backend=backend, + ) val.append(expval) mit_val = model(np.array(val).reshape(-1, 1), *optimal_params[0])[0] - # Return data - if full_output == True: + if full_output is True: return mit_val, val, optimal_params, train_val return mit_val -def calibration_matrix(nqubits, noise_model=None, nshots: int = 1000, backend=None): - """Computes the calibration matrix for readout mitigation. +def iterative_bayesian_unfolding(probabilities, response_matrix, iterations=10): + """ + Iterative Bayesian Unfolding (IBU) method for readout mitigation. + + Args: + probabilities (numpy.ndarray): the input probabilities to be unfolded. + response_matrix (numpy.ndarray): the response matrix. + iterations (int, optional): the number of iterations to perform. Defaults to 10. + + Returns: + numpy.ndarray: the unfolded probabilities. + + Reference: + 1. B. Nachman, M. Urbanek et al, *Unfolding Quantum Computer Readout Noise*. + `arXiv:1910.01969 [quant-ph] `_. + 2. S. Srinivasan, B. Pokharel et al, *Scalable Measurement Error Mitigation via Iterative Bayesian Unfolding*. + `arXiv:2210.12284 [quant-ph] `_. + """ + unfolded_probabilities = np.ones((len(probabilities), 1)) / len(probabilities) + + for _ in range(iterations): + unfolded_probabilities = unfolded_probabilities * ( + np.transpose(response_matrix) + @ (probabilities / (response_matrix @ unfolded_probabilities)) + ) + + return unfolded_probabilities + + +def get_response_matrix( + nqubits, qubit_map=None, noise_model=None, nshots: int = 10000, backend=None +): + """Computes the response matrix for readout mitigation. Args: nqubits (int): Total number of qubits. + qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``. noise_model (:class:`qibo.noise.NoiseModel`, optional): noise model used for simulating noisy computation. This matrix can be used to mitigate the effect of `qibo.noise.ReadoutError`. - nshots (int, optional): number of shots. - backend (:class:`qibo.backends.abstract.Backend`, optional): calculation engine. + nshots (int, optional): number of shots. Defaults to :math:`10000`. + backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used + in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. Returns: - numpy.ndarray : The computed (`nqubits`, `nqubits`) calibration matrix for + numpy.ndarray : the computed (`nqubits`, `nqubits`) response matrix for readout mitigation. """ - from qibo import Circuit # pylint: disable=import-outside-toplevel if backend is None: # pragma: no cover backend = GlobalBackend() - matrix = np.zeros((2**nqubits, 2**nqubits)) + response_matrix = np.zeros((2**nqubits, 2**nqubits)) for i in range(2**nqubits): - state = format(i, f"0{nqubits}b") + binary_state = format(i, f"0{nqubits}b") circuit = Circuit(nqubits, density_matrix=True) - for q, bit in enumerate(state): + for qubit, bit in enumerate(binary_state): if bit == "1": - circuit.add(gates.X(q)) + circuit.add(gates.X(qubit)) circuit.add(gates.M(*range(nqubits))) - if noise_model is not None and backend.name != "qibolab": - circuit = noise_model.apply(circuit) + circuit_result = _execute_circuit( + circuit, qubit_map, noise_model, nshots, backend=backend + ) - freq = backend.execute_circuit(circuit, nshots=nshots).frequencies() + frequencies = circuit_result.frequencies() column = np.zeros(2**nqubits) - for key in freq.keys(): - f = freq[key] / nshots - column[int(key, 2)] = f - matrix[:, i] = column + for key, value in frequencies.items(): + column[int(key, 2)] = value / nshots + response_matrix[:, i] = column - return np.linalg.inv(matrix) + return response_matrix -def apply_readout_mitigation(state, calibration_matrix): - """Updates the frequencies of the input state with the mitigated ones obtained with - ``calibration_matrix * state.frequencies()``. +def apply_resp_mat_readout_mitigation(state, response_matrix, iterations=None): + """ + Applies readout error mitigation to the given state using the provided response matrix. Args: - state (:class:`qibo.measurements.CircuitResult`): input state to be updated. - calibration_matrix (numpy.ndarray): calibration matrix for readout mitigation. + state (:class:`qibo.measurements.CircuitResult`): the input state to be updated. This state should contain the + frequencies that need to be mitigated. + response_matrix (numpy.ndarray): the response matrix for readout mitigation. + iterations (int, optional): the number of iterations to use for the Iterative Bayesian Unfolding method. + If ``None`` the 'inverse' method is used. Defaults to ``None``. Returns: - :class:`qibo.measurements.CircuitResult`: the input state with the updated frequencies. + :class:`qibo.measurements.CircuitResult`: the input state with the updated (mitigated) frequencies. """ - freq = np.zeros(2**state.nqubits) - for k, v in state.frequencies().items(): - freq[int(k, 2)] = v + frequencies = np.zeros(2 ** len(state.measurements[0].qubits)) + for key, value in state.frequencies().items(): + frequencies[int(key, 2)] = value - freq = freq.reshape(-1, 1) + frequencies = frequencies.reshape(-1, 1) - for i, val in enumerate(calibration_matrix @ freq): - state._frequencies[i] = float(val) + if iterations is None: + calibration_matrix = np.linalg.inv(response_matrix) + mitigated_frequencies = calibration_matrix @ frequencies + else: + mitigated_probabilities = iterative_bayesian_unfolding( + frequencies / np.sum(frequencies), response_matrix, iterations + ) + mitigated_frequencies = np.round( + mitigated_probabilities * np.sum(frequencies), 0 + ) + mitigated_frequencies = ( + mitigated_frequencies / np.sum(mitigated_frequencies) + ) * np.sum(frequencies) + + for i, value in enumerate(mitigated_frequencies): + state._frequencies[i] = float(value) return state def apply_randomized_readout_mitigation( - circuit, noise_model=None, nshots: int = int(1e3), ncircuits: int = 10, backend=None + circuit, + noise_model=None, + nshots: int = 10000, + ncircuits: int = 10, + qubit_map=None, + backend=None, ): - """Implements the readout mitigation method proposed in https://arxiv.org/abs/2012.09738. + """Readout mitigation method that transforms the bias in an expectation value into a measurable multiplicative factor. This factor can be eliminated at the expense of increased sampling complexity for the observable. Args: circuit (:class:`qibo.models.Circuit`): input circuit. noise_model(:class:`qibo.noise.NoiseModel`, optional): noise model used for - simulating noisy computation. This matrix can be used to mitigate the - effects of :class:`qibo.noise.ReadoutError`. - nshots (int, optional): number of shots. + simulating noisy computation. Defaults to ``None``. + nshots (int, optional): number of shots. Defaults to :math:`10000`. ncircuits (int, optional): number of randomized circuits. Each of them uses - ``int(nshots / ncircuits)`` shots. - backend (:class:`qibo.backends.abstract.Backend`): calculation engine. + ``int(nshots / ncircuits)`` shots. Defaults to 10. + qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``. + backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used + in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. Return: :class:`qibo.measurements.CircuitResult`: the state of the input circuit with mitigated frequencies. + + Reference: + 1. Ewout van den Berg, Zlatko K. Minev et al, *Model-free readout-error mitigation for quantum expectation values*. + `arXiv:2012.09738 [quant-ph] `_. """ from qibo import Circuit # pylint: disable=import-outside-toplevel from qibo.quantum_info import ( # pylint: disable=import-outside-toplevel @@ -560,7 +600,7 @@ def apply_randomized_readout_mitigation( if backend is None: # pragma: no cover backend = GlobalBackend() - qubits = circuit.queue[-1].qubits + meas_qubits = circuit.measurements[0].qubits nshots_r = int(nshots / ncircuits) freq = np.zeros((ncircuits, 2), object) for k in range(ncircuits): @@ -568,11 +608,11 @@ def apply_randomized_readout_mitigation( circuit_c.queue.pop() cal_circuit = Circuit(circuit.nqubits, density_matrix=True) - x_gate = random_pauli(len(qubits), 1, subset=["I", "X"]).queue + x_gate = random_pauli(circuit.nqubits, 1, subset=["I", "X"]).queue error_map = {} - for gate in x_gate: - if gate.name == "x": + for j, gate in enumerate(x_gate): + if isinstance(gate, gates.X) and gate.qubits[0] in meas_qubits: error_map[gate.qubits[0]] = 1 circuits = [circuit_c, cal_circuit] @@ -580,10 +620,11 @@ def apply_randomized_readout_mitigation( freqs = [] for circ in circuits: circ.add(x_gate) - circ.add(gates.M(*qubits)) - if noise_model is not None and backend.name != "qibolab": - circ = noise_model.apply(circ) - result = backend.execute_circuit(circ, nshots=nshots_r) + circ.add(gates.M(*meas_qubits)) + + result = _execute_circuit( + circ, qubit_map, noise_model, nshots_r, backend=backend + ) result._samples = result.apply_bitflips(error_map) results.append(result) freqs.append(result.frequencies(binary=False)) @@ -592,8 +633,345 @@ def apply_randomized_readout_mitigation( for j in range(2): results[j].nshots = nshots freq_sum = freq[0, j] - for f in freq[1::, j]: - freq_sum += f + for frs in freq[1::, j]: + freq_sum += frs results[j]._frequencies = freq_sum return results + + +def get_expectation_val_with_readout_mitigation( + circuit, + observable, + noise_model=None, + nshots: int = 10000, + readout=None, + qubit_map=None, + backend=None, +): + """ + Applies readout error mitigation to the given circuit and observable. + + Args: + circuit (qibo.models.Circuit): input circuit. + observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): The observable to be measured. + noise_model (qibo.models.noise.Noise, optional): the noise model to be applied. Defaults to ``None``. + nshots (int, optional): the number of shots for the circuit execution. Defaults to :math:`10000`. + readout (dict, optional): a dictionary that may contain the following keys: + + * ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation. + * response_matrix: numpy.ndarray, used for applying a pre-computed response matrix for readout error mitigation. + * ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}. + + qubit_map (list, optional): the qubit map. If None, a list of range of circuit's qubits is used. Defaults to ``None``. + backend (qibo.backends.abstract.Backend, optional): the backend to be used in the execution. + If None, it uses the global backend. Defaults to ``None``. + + Returns: + float: the mitigated expectation value of the observable. + """ + if backend is None: # pragma: no cover + backend = GlobalBackend() + if readout is None: # pragma: no cover + readout = {} + + if "ncircuits" in readout: + circuit_result, circuit_result_cal = apply_randomized_readout_mitigation( + circuit, noise_model, nshots, readout["ncircuits"], backend + ) + else: + circuit_result = _execute_circuit( + circuit, qubit_map, noise_model, nshots, backend=backend + ) + if "response_matrix" in readout: + circuit_result = apply_resp_mat_readout_mitigation( + circuit_result, + readout["response_matrix"], + readout.get("ibu_iters", None), + ) + + exp_val = circuit_result.expectation_from_samples(observable) + + if "ncircuits" in readout: + exp_val /= circuit_result_cal.expectation_from_samples(observable) + + return exp_val + + +def sample_clifford_training_circuit( + circuit, + backend=None, +): + """Samples a training circuit for CDR by susbtituting all the non-Clifford gates. + + Args: + circuit (:class:`qibo.models.Circuit`): circuit to sample from. + backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used + in the execution. If ``None``, it uses :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. + + Returns: + :class:`qibo.models.Circuit`: the sampled circuit. + """ + from qibo.quantum_info import ( # pylint: disable=import-outside-toplevel + random_clifford, + ) + + if backend is None: # pragma: no cover + backend = GlobalBackend() + + non_clifford_gates_indices = [ + i + for i, gate in enumerate(circuit.queue) + if not gate.clifford and not isinstance(gate, gates.M) + ] + + if not non_clifford_gates_indices: + raise_error(ValueError, "No non-Clifford gate found, no circuit sampled.") + + sampled_circuit = circuit.__class__(**circuit.init_kwargs) + + for i, gate in enumerate(circuit.queue): + if isinstance(gate, gates.M): + for q in gate.qubits: + gate_rand = gates.Unitary( + random_clifford(1, backend=backend, return_circuit=False), + q, + ) + gate_rand.clifford = True + sampled_circuit.add(gate_rand) + sampled_circuit.add(gate) + else: + if i in non_clifford_gates_indices: + gate = gates.Unitary( + random_clifford( + len(gate.qubits), backend=backend, return_circuit=False + ), + *gate.qubits, + ) + gate.clifford = True + sampled_circuit.add(gate) + + return sampled_circuit + + +def error_sensitive_circuit(circuit, observable, backend=None): + """ + Generates a Clifford circuit that preserves the same circuit frame as the input circuit, and stabilizes the specified Pauli observable. + + Args: + circuit (:class:`qibo.models.Circuit`): input circuit. + observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): Pauli observable to be measured. + backend (:class:`qibo.backends.abstract.Backend`, optional): backend to be used + in the execution. if ``None``, it uses :class:`qibo.backends.GlobalBackend`. + Defaults to ``None``. + + Returns: + :class:`qibo.models.Circuit`: the error sensitive circuit. + :class:`qibo.models.Circuit`: the sampled Clifford circuit. + list: the list of adjustment gates. + + Reference: + 1. Dayue Qin, Yanzhu Chen et al, *Error statistics and scalability of quantum error mitigation formulas*. + `arXiv:2112.06255 [quant-ph] `_. + """ + from qibo.quantum_info import ( # pylint: disable=import-outside-toplevel + comp_basis_to_pauli, + random_clifford, + vectorization, + ) + + if backend is None: # pragma: no cover + backend = GlobalBackend() + + sampled_circuit = sample_clifford_training_circuit(circuit, backend=backend) + unitary_matrix = sampled_circuit.unitary(backend=backend) + num_qubits = sampled_circuit.nqubits + + comp_to_pauli = comp_basis_to_pauli(num_qubits, backend=backend) + observable.nqubits = num_qubits + observable_liouville = vectorization( + np.transpose(np.conjugate(unitary_matrix)) @ observable.matrix @ unitary_matrix, + order="row", + backend=backend, + ) + observable_pauli_liouville = comp_to_pauli @ observable_liouville + + index = int(np.where(abs(observable_pauli_liouville) >= 1e-5)[0][0]) + + observable_pauli = list(product(["I", "X", "Y", "Z"], repeat=num_qubits))[index] + + pauli_gates = { + "I": backend.cast(gates.I(0).matrix(backend=backend)), + "X": backend.cast(gates.X(0).matrix(backend=backend)), + "Y": backend.cast(gates.Y(0).matrix(backend=backend)), + "Z": backend.cast(gates.Z(0).matrix(backend=backend)), + } + + adjustment_gates = [] + for i in range(num_qubits): + observable_i = pauli_gates[observable_pauli[i]] + random_init = pauli_gates["I"] + while np.any(abs(observable_i - pauli_gates["Z"]) > 1e-5) and np.any( + abs(observable_i - pauli_gates["I"]) > 1e-5 + ): + random_init = random_clifford(1, backend=backend, return_circuit=False) + observable_i = ( + np.conjugate(np.transpose(random_init)) + @ pauli_gates[observable_pauli[i]] + @ random_init + ) + + adjustment_gate = gates.Unitary(random_init, i) + adjustment_gate.clifford = True + adjustment_gates.append(adjustment_gate) + + sensitive_circuit = sampled_circuit.__class__(**sampled_circuit.init_kwargs) + + for gate in adjustment_gates: + sensitive_circuit.add(gate) + for gate in sampled_circuit.queue: + sensitive_circuit.add(gate) + + return sensitive_circuit, sampled_circuit, adjustment_gates + + +def ICS( + circuit, + observable, + readout=None, + qubit_map=None, + noise_model=None, + nshots=int(1e4), + n_training_samples=10, + full_output=False, + backend=None, +): + """ + Computes the Important Clifford Sampling method. + + Args: + circuit (:class:`qibo.models.Circuit`): input circuit. + observable (:class:`qibo.hamiltonians.Hamiltonian/:class:`qibo.hamiltonians.SymbolicHamiltonian`): the observable to be measured. + readout (dict, optional): a dictionary that may contain the following keys: + + * ncircuits: int, specifies the number of random circuits to use for the randomized method of readout error mitigation. + * response_matrix: numpy.ndarray, used for applying a pre-computed response matrix for readout error mitigation. + * ibu_iters: int, specifies the number of iterations for the iterative Bayesian unfolding method of readout error mitigation. If provided, the corresponding readout error mitigation method is used. Defaults to {}. + + qubit_map (list, optional): the qubit map. If ``None``, a list of range of circuit's qubits is used. Defaults to ``None``. + noise_model (qibo.models.noise.Noise, optional): the noise model to be applied. Defaults to ``None``. + nshots (int, optional): the number of shots for the circuit execution. Defaults to :math:`10000`. + n_training_samples (int, optional): the number of training samples. Defaults to 10. + full_output (bool, optional): if ``True``, this function returns additional + information: ``val``, ``optimal_params``, ``train_val``. Defaults to ``False``. + backend (qibo.backends.abstract.Backend, optional): the backend to be used in the execution. + If None, it uses the global backend. Defaults to ``None``. + + Returns: + mitigated_expectation (float): the mitigated expectated value. + mitigated_expectation_std (float): the standard deviation of the mitigated expectated value. + dep_param (float): the depolarizing parameter. + dep_param_std (float): the standard deviation of the depolarizing parameter. + lambda_list (list): the list of the depolarizing parameters. + data (dict): the data dictionary containing the noise-free and noisy expectation values obtained with the training circuits. + + Reference: + 1. Dayue Qin, Yanzhu Chen et al, *Error statistics and scalability of quantum error mitigation formulas*. + `arXiv:2112.06255 [quant-ph] `_. + """ + if backend is None: # pragma: no cover + backend = GlobalBackend() + if readout is None: + readout = {} + + if qubit_map is None: + qubit_map = list(range(circuit.nqubits)) + + training_circuits = [ + error_sensitive_circuit(circuit, observable, backend=backend)[0] + for _ in range(n_training_samples) + ] + + data = {"noise-free": [], "noisy": []} + lambda_list = [] + + for training_circuit in training_circuits: + circuit_result = training_circuit(nshots=nshots) + expectation = observable.expectation_from_samples(circuit_result.frequencies()) + + noisy_expectation = get_expectation_val_with_readout_mitigation( + training_circuit, + observable, + noise_model, + nshots, + readout, + qubit_map, + backend=backend, + ) + + data["noise-free"].append(expectation) + data["noisy"].append(noisy_expectation) + lambda_list.append(1 - noisy_expectation / expectation) + + dep_param = np.mean(lambda_list) + dep_param_std = np.std(lambda_list) + + noisy_expectation = get_expectation_val_with_readout_mitigation( + circuit, observable, noise_model, nshots, readout, qubit_map, backend=backend + ) + one_dep_squared = (1 - dep_param) ** 2 + dep_std_squared = dep_param_std**2 + + mitigated_expectation = ( + (1 - dep_param) * noisy_expectation / (one_dep_squared + dep_std_squared) + ) + mitigated_expectation_std = ( + dep_param_std + * abs(noisy_expectation) + * abs((1 - dep_param) ** 2 - dep_std_squared) + / (one_dep_squared + dep_std_squared) ** 2 + ) + + if full_output: + return ( + mitigated_expectation, + mitigated_expectation_std, + dep_param, + dep_param_std, + lambda_list, + data, + ) + + return mitigated_expectation + + +def _execute_circuit(circuit, qubit_map, noise_model=None, nshots=10000, backend=None): + """ + Helper function to execute the given circuit with the specified parameters. + + Args: + circuit (qibo.models.Circuit): input circuit. + qubit_map (list): the qubit map. If ``None``, a list of range of circuit's qubits is used. Defaults to ``None``. + noise_model (qibo.models.noise.Noise, optional): The noise model to be applied. Defaults to ``None``. + nshots (int): the number of shots for the circuit execution. Defaults to :math:`10000`.. + backend (qibo.backends.abstract.Backend, optional): the backend to be used in the execution. + If None, it uses the global backend. Defaults to ``None``. + + Returns: + qibo.states.CircuitResult: The result of the circuit execution. + """ + from qibo.transpiler.placer import Custom + + if backend is None: # pragma: no cover + backend = GlobalBackend() + elif backend.name == "qibolab": # pragma: no cover + backend.transpiler.passes[1] = Custom( + map=qubit_map, connectivity=backend.platform.topology + ) + elif noise_model is not None: + circuit = noise_model.apply(circuit) + + circuit_result = backend.execute_circuit(circuit, nshots=nshots) + + return circuit_result diff --git a/tests/test_models_error_mitigation.py b/tests/test_models_error_mitigation.py index d602d0ae4c..035daabd1b 100644 --- a/tests/test_models_error_mitigation.py +++ b/tests/test_models_error_mitigation.py @@ -2,14 +2,16 @@ import pytest from qibo import Circuit, gates +from qibo.backends import construct_backend from qibo.hamiltonians import SymbolicHamiltonian from qibo.models.error_mitigation import ( CDR, + ICS, ZNE, - apply_randomized_readout_mitigation, - apply_readout_mitigation, - calibration_matrix, - sample_training_circuit, + get_expectation_val_with_readout_mitigation, + get_response_matrix, + sample_clifford_training_circuit, + sample_training_circuit_cdr, vnCDR, ) from qibo.noise import DepolarizingError, NoiseModel, ReadoutError @@ -17,15 +19,18 @@ from qibo.symbols import Z -def get_noise_model(error, gate, cal_matrix=[False, None]): +def get_noise_model(error, gate, resp_matrix=[False, None]): noise = NoiseModel() noise.add(error, gate) - if cal_matrix[0]: - noise.add(ReadoutError(probabilities=cal_matrix[1]), gate=gates.M) + if resp_matrix[0]: + noise.add(ReadoutError(probabilities=resp_matrix[1]), gate=gates.M) + return noise -def get_circuit(nqubits): +def get_circuit(nqubits, nmeas=None): + if nmeas is None: + nmeas = nqubits # Define the circuit hz = 0.5 hx = 0.5 @@ -42,31 +47,30 @@ def get_circuit(nqubits): c.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2)) c.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(1, nqubits, 2)) c.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2)) - c.add(gates.M(*range(nqubits))) - return c + c.add(gates.M(*range(nmeas))) + return c -from qibo.backends import construct_backend backend = construct_backend("numpy") -# # Generate random calibration matrices -cal_matrix_1q = random_stochastic_matrix( +# # Generate random response matrices +resp_matrix_1q = random_stochastic_matrix( 2, diagonally_dominant=True, seed=2, backend=backend ) -cal_matrix_3q = random_stochastic_matrix( - 8, diagonally_dominant=True, seed=2, backend=backend +resp_matrix_2q = random_stochastic_matrix( + 4, diagonally_dominant=True, seed=2, backend=backend ) @pytest.mark.parametrize( "nqubits,noise,insertion_gate,readout", [ - (3, get_noise_model(DepolarizingError(0.1), gates.CNOT), "CNOT", {}), + (3, get_noise_model(DepolarizingError(0.1), gates.CNOT), "CNOT", None), ( 3, get_noise_model(DepolarizingError(0.1), gates.CNOT), "CNOT", - {"calibration_matrix": cal_matrix_3q}, + {"response_matrix": resp_matrix_2q, "ibu_iters": None}, ), ( 3, @@ -74,12 +78,18 @@ def get_circuit(nqubits): "CNOT", {"ncircuits": 2}, ), - (1, get_noise_model(DepolarizingError(0.1), gates.RX), "RX", {}), + (1, get_noise_model(DepolarizingError(0.1), gates.RX), "RX", None), + ( + 1, + get_noise_model(DepolarizingError(0.3), gates.RX), + "RX", + {"response_matrix": resp_matrix_1q, "ibu_iters": None}, + ), ( 1, get_noise_model(DepolarizingError(0.3), gates.RX), "RX", - {"calibration_matrix": cal_matrix_1q}, + {"response_matrix": resp_matrix_1q, "ibu_iters": 10}, ), (1, get_noise_model(DepolarizingError(0.1), gates.RX), "RX", {"ncircuits": 2}), ], @@ -94,13 +104,19 @@ def test_zne(backend, nqubits, noise, solve, insertion_gate, readout): tf.config.threading.set_intra_op_parallelism_threads = 1 else: backend.set_threads(1) + + if nqubits == 1: + nmeas = 1 + else: + nmeas = nqubits - 1 # Define the circuit - c = get_circuit(nqubits) + c = get_circuit(nqubits, nmeas) # Define the observable - obs = np.prod([Z(i) for i in range(nqubits)]) + obs = np.prod([Z(i) for i in range(nmeas)]) + obs_exact = SymbolicHamiltonian(obs, nqubits=nqubits, backend=backend) obs = SymbolicHamiltonian(obs, backend=backend) # Noise-free expected value - exact = obs.expectation(backend.execute_circuit(c).state()) + exact = obs_exact.expectation(backend.execute_circuit(c).state()) # Noisy expected value without mitigation state = backend.execute_circuit(noise.apply(c), nshots=10000) noisy = state.expectation_from_samples(obs) @@ -116,6 +132,7 @@ def test_zne(backend, nqubits, noise, solve, insertion_gate, readout): readout=readout, backend=backend, ) + assert np.abs(exact - estimate) <= np.abs(exact - noisy) @@ -124,18 +141,23 @@ def test_zne(backend, nqubits, noise, solve, insertion_gate, readout): @pytest.mark.parametrize( "noise,readout", [ - (get_noise_model(DepolarizingError(0.1), gates.CNOT), {}), + (get_noise_model(DepolarizingError(0.1), gates.CNOT), None), + ( + get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]), + {"response_matrix": resp_matrix_2q, "ibu_iters": None}, + ), ( - get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]), - {"calibration_matrix": cal_matrix_3q}, + get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]), + {"response_matrix": resp_matrix_2q, "ibu_iters": 10}, ), ( - get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]), + get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]), {"ncircuits": 2}, ), ], ) def test_cdr(backend, nqubits, noise, full_output, readout): + """Test that CDR reduces the noise.""" if backend.name == "tensorflow": import tensorflow as tf @@ -143,14 +165,19 @@ def test_cdr(backend, nqubits, noise, full_output, readout): tf.config.threading.set_intra_op_parallelism_threads = 1 else: backend.set_threads(1) - """Test that CDR reduces the noise.""" + + if nqubits == 1: + nmeas = 1 + else: + nmeas = nqubits - 1 # Define the circuit - c = get_circuit(nqubits) + c = get_circuit(nqubits, nmeas) # Define the observable - obs = np.prod([Z(i) for i in range(nqubits)]) + obs = np.prod([Z(i) for i in range(nmeas)]) + obs_exact = SymbolicHamiltonian(obs, nqubits=nqubits, backend=backend) obs = SymbolicHamiltonian(obs, backend=backend) # Noise-free expected value - exact = obs.expectation(backend.execute_circuit(c).state()) + exact = obs_exact.expectation(backend.execute_circuit(c).state()) # Noisy expected value without mitigation state = backend.execute_circuit(noise.apply(c), nshots=10000) noisy = state.expectation_from_samples(obs) @@ -167,6 +194,7 @@ def test_cdr(backend, nqubits, noise, full_output, readout): ) if full_output: estimate = estimate[0] + assert np.abs(exact - estimate) <= np.abs(exact - noisy) @@ -189,36 +217,28 @@ def test_sample_training_circuit(nqubits): c.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(1, nqubits, 2)) c.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2)) c.add(gates.M(q) for q in range(nqubits)) + + for j in range(len(c.queue)): + c.queue[j].clifford = True + with pytest.raises(ValueError): + sample_training_circuit_cdr(c) with pytest.raises(ValueError): - sample_training_circuit(c) + sample_clifford_training_circuit(c) @pytest.mark.parametrize( "nqubits,noise,insertion_gate,readout", [ - (3, get_noise_model(DepolarizingError(0.1), gates.CNOT), "CNOT", {}), - ( - 3, - get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]), - "CNOT", - {"calibration_matrix": cal_matrix_3q}, - ), - ( - 3, - get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, cal_matrix_3q]), - "CNOT", - {"ncircuits": 2}, - ), - (1, get_noise_model(DepolarizingError(0.1), gates.RX), "RX", {}), + (1, get_noise_model(DepolarizingError(0.1), gates.RX), "RX", None), ( 1, - get_noise_model(DepolarizingError(0.1), gates.RX, [True, cal_matrix_1q]), + get_noise_model(DepolarizingError(0.1), gates.RX, [True, resp_matrix_1q]), "RX", - {"calibration_matrix": cal_matrix_1q}, + {"response_matrix": resp_matrix_1q, "ibu_iters": 10}, ), ( 1, - get_noise_model(DepolarizingError(0.1), gates.RX, [True, cal_matrix_1q]), + get_noise_model(DepolarizingError(0.1), gates.RX, [True, resp_matrix_1q]), "RX", {"ncircuits": 2}, ), @@ -242,12 +262,6 @@ def test_vncdr(backend, nqubits, noise, full_output, insertion_gate, readout): # Noise-free expected value exact = obs.expectation(backend.execute_circuit(c).state()) # Noisy expected value without mitigation - if "calibration_matrix" in readout.keys() or "ncircuits" in readout.keys(): - if nqubits == 1: - p = cal_matrix_1q - elif nqubits == 3: - p = cal_matrix_3q - # noise.add(ReadoutError(probabilities=p),gate=gates.M) state = backend.execute_circuit(noise.apply(c), nshots=10000) noisy = state.expectation_from_samples(obs) # Mitigated expected value @@ -265,12 +279,14 @@ def test_vncdr(backend, nqubits, noise, full_output, insertion_gate, readout): ) if full_output: estimate = estimate[0] + assert np.abs(exact - estimate) <= np.abs(exact - noisy) -@pytest.mark.parametrize("nqubits", [3]) -@pytest.mark.parametrize("method", ["cal_matrix", "temme"]) -def test_readout_mitigation(backend, nqubits, method): +@pytest.mark.parametrize("nqubits,nmeas", [(3, 2)]) +@pytest.mark.parametrize("method", ["response_matrix", "randomized"]) +@pytest.mark.parametrize("ibu_iters", [None, 10]) +def test_readout_mitigation(backend, nqubits, nmeas, method, ibu_iters): if backend.name == "tensorflow": import tensorflow as tf @@ -279,30 +295,86 @@ def test_readout_mitigation(backend, nqubits, method): else: backend.set_threads(1) nshots = 10000 - p = random_stochastic_matrix(2**nqubits, diagonally_dominant=True, seed=5) + p = random_stochastic_matrix(2**nmeas, diagonally_dominant=True, seed=5) noise = NoiseModel() noise.add(ReadoutError(probabilities=p), gate=gates.M) - if method == "cal_matrix": - calibration = calibration_matrix(nqubits, noise, nshots=nshots, backend=backend) + if method == "response_matrix": + response = get_response_matrix( + nmeas, None, noise, nshots=nshots, backend=backend + ) + readout = {"response_matrix": response, "ibu_iters": ibu_iters} + elif method == "randomized": + readout = {"ncircuits": 10} # Define the observable - obs = np.prod([Z(i) for i in range(nqubits)]) + obs = np.prod([Z(i) for i in range(nmeas)]) obs = SymbolicHamiltonian(obs, backend=backend) # get noise free expected val - c = get_circuit(nqubits) + c = get_circuit(nqubits, nmeas) true_state = backend.execute_circuit(c, nshots=nshots) true_val = true_state.expectation_from_samples(obs) # get noisy expected val state = backend.execute_circuit(noise.apply(c), nshots=nshots) noisy_val = state.expectation_from_samples(obs) - if method == "cal_matrix": - mit_state = apply_readout_mitigation(state, calibration) - mit_val = mit_state.expectation_from_samples(obs) - elif method == "temme": - ncircuits = 10 - result, result_cal = apply_randomized_readout_mitigation( - c, noise, nshots, ncircuits, backend - ) - mit_val = result.expectation_from_samples( - obs - ) / result_cal.expectation_from_samples(obs) + + mit_val = get_expectation_val_with_readout_mitigation( + c, obs, noise, nshots, readout, backend=backend + ) + assert np.abs(true_val - mit_val) <= np.abs(true_val - noisy_val) + + +@pytest.mark.parametrize("nqubits", [3]) +@pytest.mark.parametrize("full_output", [False, True]) +@pytest.mark.parametrize( + "noise,readout", + [ + (get_noise_model(DepolarizingError(0.1), gates.CNOT), None), + ( + get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]), + {"response_matrix": resp_matrix_2q, "ibu_iters": None}, + ), + ( + get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]), + {"response_matrix": resp_matrix_2q, "ibu_iters": 10}, + ), + ( + get_noise_model(DepolarizingError(0.1), gates.CNOT, [True, resp_matrix_2q]), + {"ncircuits": 2}, + ), + ], +) +def test_ics(backend, nqubits, noise, full_output, readout): + if backend.name == "tensorflow": + import tensorflow as tf + + tf.config.threading.set_inter_op_parallelism_threads = 1 + tf.config.threading.set_intra_op_parallelism_threads = 1 + else: + backend.set_threads(1) + """Test that ICS reduces the noise.""" + # Define the circuit + c = get_circuit(nqubits, nqubits - 1) + # Define the observable + obs = np.prod([Z(i) for i in range(nqubits - 1)]) + obs_exact = SymbolicHamiltonian(obs, nqubits=nqubits, backend=backend) + obs = SymbolicHamiltonian(obs, backend=backend) + # Noise-free expected value + exact = obs_exact.expectation(backend.execute_circuit(c).state()) + # Noisy expected value without mitigation + state = backend.execute_circuit(noise.apply(c), nshots=10000) + noisy = state.expectation_from_samples(obs) + # Mitigated expected value + estimate = ICS( + circuit=c, + observable=obs, + noise_model=noise, + nshots=10000, + n_training_samples=20, + full_output=full_output, + readout=readout, + backend=backend, + ) + if full_output: + estimate = estimate[0] + + assert np.abs(exact - estimate) <= np.abs(exact - noisy)