Merge pull request #1129 from qiboteam/cdr

Quantum Error Mitigation: Readout, Importance Clifford Sampling

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 <https://arxiv.org/abs/1910.01969>`_ 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 <https://arxiv.org/abs/2103.12680>`_ 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 <https://arxiv.org/abs/2103.12680>`_ 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 <https://arxiv.org/abs/2103.12680>`_ 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

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 <error-mitigation>`.
 
 .. _timeevol-example:

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