replace baseclass of IRB analysis.

qiskit_experiments/curve_analysis/__init__.py

@@ -90,6 +90,7 @@
     guess.constant_sinusoidal_offset
     guess.constant_spectral_offset
     guess.exp_decay
+    guess.rb_decay
     guess.full_width_half_max
     guess.frequency
     guess.max_height

qiskit_experiments/curve_analysis/guess.py

@@ -345,3 +345,67 @@ def constant_sinusoidal_offset(y: np.ndarray) -> float:
     minv, _ = min_height(y, percentile=5)
 
     return 0.5 * (maxv + minv)
+
+
+def rb_decay(
+    x: np.ndarray,
+    y: np.ndarray,
+    b: float = 0.5,
+    a: Optional[float] = None,
+) -> float:
+    r"""Get base of exponential decay function which is assumed to be close to 1.
+
+    This assumes following model:
+
+    .. math::
+
+        y(x) = a \alpha^x + b.
+
+    This model can be often seen in randomized benchmarking analysis.
+    Now we assume :math:`\alpha \simeq 1` and replace it with :math:`\alpha = 1 - p`
+    where :math:`p \ll 1`. Then, above function can be approximated to
+
+    .. math::
+
+        y(x) = a \left(1 - px + \frac{x(x-1)}{2} p^2 + ...\right) + b,
+
+    here we assume :math:`px \ll 1` to ignore higher order terms of :math:`p`.
+    Likely parameter :math:`p` of interest can be solved as
+
+    .. math::
+
+        p(x) = \frac{a + b - y(x)}{ax}.
+
+    The :math:`p(x)` is averaged with the weight of :math:`1/x`.
+    When estimated :math:`\alpha = 1-p < 0.9`, this calcualtes following value instead.
+
+    .. math::
+
+        \alpha \simeq dy ^ {1 / dx},
+
+    where :math:`dx = x[1] - x[0]` and :math:`dy = (y[1] - b) / (y[0] - b)`.
+
+    Args:
+        x: Array of x values.
+        y: Array of y values.
+        a: Coefficient of decay function. If not provided this is detemined so that
+            :math:`y(0) = 1.0`.
+        b: Asymptote of decay function.
+
+    Returns:
+         Base of decay function.
+    """
+    if a is None:
+        a = 1 - b
+
+    a_guess_tmp = np.average(1 - (a + b - y) / (a * x), weights=1 / x)
+
+    if a_guess_tmp > 0.9:
+        return a_guess_tmp
+
+    # Violate assumption, then switch to the conventional method
+    # Use the first two points to guess the decay param
+    dx = x[1] - x[0]
+    dy = (y[1] - b) / (y[0] - b)
+
+    return dy ** (1 / dx)

qiskit_experiments/library/randomized_benchmarking/interleaved_rb_analysis.py

@@ -18,10 +18,8 @@
 import qiskit_experiments.curve_analysis as curve
 from qiskit_experiments.framework import AnalysisResultData
 
-from .rb_analysis import RBAnalysis
 
-
-class InterleavedRBAnalysis(RBAnalysis):
+class InterleavedRBAnalysis(curve.CurveAnalysis):
     r"""A class to analyze interleaved randomized benchmarking experiment.
 
     # section: overview
@@ -69,22 +67,20 @@ class InterleavedRBAnalysis(RBAnalysis):
     # section: fit_parameters
         defpar a:
             desc: Height of decay curve.
-            init_guess: Determined by the average :math:`a` of the standard and interleaved RB.
+            init_guess: Determined by :math:`1 - b`.
             bounds: [0, 1]
         defpar b:
             desc: Base line.
-            init_guess: Determined by the average :math:`b` of the standard and interleaved RB.
-                Usually equivalent to :math:`(1/2)^n` where :math:`n` is number of qubit.
+            init_guess: Determined by :math:`(1/2)^n` where :math:`n` is number of qubit.
             bounds: [0, 1]
         defpar \alpha:
             desc: Depolarizing parameter.
-            init_guess: Determined by the slope of :math:`(y - b)^{-x}` of the first and the
-                second data point of the standard RB.
+            init_guess: Determined by :func:`~rb_decay` with standard RB curve.
             bounds: [0, 1]
         defpar \alpha_c:
             desc: Ratio of the depolarizing parameter of interleaved RB to standard RB curve.
-            init_guess: Estimate :math:`\alpha' = \alpha_c \alpha` from the
-                interleaved RB curve, then divide this by the initial guess of :math:`\alpha`.
+            init_guess: Determined by alpha of interleaved RB curve devided by one of
+                standard RB curve. Both alpha values are estimated by :func:`~rb_decay`.
             bounds: [0, 1]
 
     # section: reference
@@ -120,7 +116,6 @@ def _default_options(cls):
         """Default analysis options."""
         default_options = super()._default_options()
         default_options.result_parameters = ["alpha", "alpha_c"]
-        default_options.gate_error_ratio = False
         return default_options
 
     def _generate_fit_guesses(
@@ -141,27 +136,60 @@ def _generate_fit_guesses(
             b=(0, 1),
         )
 
+        b_guess = 1 / 2**self._num_qubits
+        a_guess = 1 - b_guess
+
         # for standard RB curve
         std_curve = self._data(series_name="Standard")
-        opt_std = user_opt.copy()
-        opt_std = self._initial_guess(opt_std, std_curve.x, std_curve.y, self._num_qubits)
+        alpha_std = curve.guess.rb_decay(std_curve.x, std_curve.y, a=a_guess, b=b_guess)
 
         # for interleaved RB curve
         int_curve = self._data(series_name="Interleaved")
-        opt_int = user_opt.copy()
-        if opt_int.p0["alpha_c"] is not None:
-            opt_int.p0["alpha"] = opt_std.p0["alpha"] * opt_int.p0["alpha_c"]
-        opt_int = self._initial_guess(opt_int, int_curve.x, int_curve.y, self._num_qubits)
+        alpha_int = curve.guess.rb_decay(int_curve.x, int_curve.y, a=a_guess, b=b_guess)
+
+        alpha_c = min(alpha_int / alpha_std, 1.0)
 
         user_opt.p0.set_if_empty(
-            a=np.mean([opt_std.p0["a"], opt_int.p0["a"]]),
-            alpha=opt_std.p0["alpha"],
-            alpha_c=min(opt_int.p0["alpha"] / opt_std.p0["alpha"], 1),
-            b=np.mean([opt_std.p0["b"], opt_int.p0["b"]]),
+            b=b_guess,
+            a=a_guess,
+            alpha=alpha_std,
+            alpha_c=alpha_c,
         )
 
         return user_opt
 
+    def _format_data(self, data: curve.CurveData) -> curve.CurveData:
+        """Data format with averaging with sampling strategy."""
+        # TODO Eventually move this to data processor, then create RB data processor.
+
+        # take average over the same x value by regenerating sigma from variance of y values
+        series, xdata, ydata, sigma, shots = curve.data_processing.multi_mean_xy_data(
+            series=data.data_index,
+            xdata=data.x,
+            ydata=data.y,
+            sigma=data.y_err,
+            shots=data.shots,
+            method="sample",
+        )
+
+        # sort by x value in ascending order
+        series, xdata, ydata, sigma, shots = curve.data_processing.data_sort(
+            series=series,
+            xdata=xdata,
+            ydata=ydata,
+            sigma=sigma,
+            shots=shots,
+        )
+
+        return curve.CurveData(
+            label="fit_ready",
+            x=xdata,
+            y=ydata,
+            y_err=sigma,
+            shots=shots,
+            data_index=series,
+        )
+
     def _extra_database_entry(self, fit_data: curve.FitData) -> List[AnalysisResultData]:
         """Calculate EPC."""
         nrb = 2**self._num_qubits

qiskit_experiments/library/randomized_benchmarking/rb_analysis.py

@@ -16,11 +16,9 @@
 from collections import defaultdict
 from typing import Dict, List, Sequence, Tuple, Union, Optional, TYPE_CHECKING
 
-import numpy as np
 from qiskit.exceptions import QiskitError
 
 import qiskit_experiments.curve_analysis as curve
-from qiskit_experiments.curve_analysis.data_processing import multi_mean_xy_data, data_sort
 from qiskit_experiments.exceptions import AnalysisError
 from qiskit_experiments.framework import AnalysisResultData, ExperimentData
 from qiskit_experiments.database_service import DbAnalysisResultV1
@@ -54,17 +52,15 @@ class RBAnalysis(curve.CurveAnalysis):
     # section: fit_parameters
         defpar a:
             desc: Height of decay curve.
-            init_guess: Determined by :math:`(y - b) / \alpha^x`.
+            init_guess: Determined by :math:`1 - b`.
             bounds: [0, 1]
         defpar b:
             desc: Base line.
-            init_guess: Determined by the average :math:`b` of the standard and interleaved RB.
-                Usually equivalent to :math:`(1/2)^n` where :math:`n` is number of qubit.
+            init_guess: Determined by :math:`(1/2)^n` where :math:`n` is number of qubit.
             bounds: [0, 1]
         defpar \alpha:
             desc: Depolarizing parameter.
-            init_guess: Determined by the slope of :math:`(y - b)^{-x}` of the first and the
-                second data point.
+            init_guess: Determined by :func:`~rb_decay`.
             bounds: [0, 1]
 
     # section: reference
@@ -139,34 +135,24 @@ def _generate_fit_guesses(
             b=(0, 1),
         )
 
-        return self._initial_guess(user_opt, curve_data.x, curve_data.y, self._num_qubits)
+        b_guess = 1 / 2**self._num_qubits
+        a_guess = 1 - b_guess
+        alpha_guess = curve.guess.rb_decay(curve_data.x, curve_data.y, a=a_guess, b=b_guess)
 
-    @staticmethod
-    def _initial_guess(
-        opt: curve.FitOptions, x_values: np.ndarray, y_values: np.ndarray, num_qubits: int
-    ) -> curve.FitOptions:
-        """Create initial guess with experiment data."""
-        opt.p0.set_if_empty(b=1 / 2**num_qubits)
-
-        # Use the first two points to guess the decay param
-        dcliff = x_values[1] - x_values[0]
-        dy = (y_values[1] - opt.p0["b"]) / (y_values[0] - opt.p0["b"])
-        alpha_guess = dy ** (1 / dcliff)
-
-        opt.p0.set_if_empty(alpha=alpha_guess if alpha_guess < 1.0 else 0.99)
-
-        if y_values[0] > opt.p0["b"]:
-            opt.p0.set_if_empty(a=(y_values[0] - opt.p0["b"]) / (opt.p0["alpha"] ** x_values[0]))
-        else:
-            opt.p0.set_if_empty(a=0.95)
+        user_opt.p0.set_if_empty(
+            b=b_guess,
+            a=a_guess,
+            alpha=alpha_guess,
+        )
 
-        return opt
+        return user_opt
 
     def _format_data(self, data: curve.CurveData) -> curve.CurveData:
         """Data format with averaging with sampling strategy."""
+        # TODO Eventually move this to data processor, then create RB data processor.
 
         # take average over the same x value by regenerating sigma from variance of y values
-        series, xdata, ydata, sigma, shots = multi_mean_xy_data(
+        series, xdata, ydata, sigma, shots = curve.data_processing.multi_mean_xy_data(
             series=data.data_index,
             xdata=data.x,
             ydata=data.y,
@@ -176,7 +162,7 @@ def _format_data(self, data: curve.CurveData) -> curve.CurveData:
         )
 
         # sort by x value in ascending order
-        series, xdata, ydata, sigma, shots = data_sort(
+        series, xdata, ydata, sigma, shots = curve.data_processing.data_sort(
             series=series,
             xdata=xdata,
             ydata=ydata,

releasenotes/notes/cleanup-rb-experiment-f17b6e674ae4e473.yaml

@@ -1,4 +1,8 @@
 ---
+features:
+  - |
+    Curve fit parameter guess function :func:`~rb_decay` has been added. 
+    This improves inital parameter estimation of randomized benchmark experiments.
 upgrade:
   - |
     Computation of error per gates (EPGs) from EPC in :class:`RBAnalysis` has been upgraded.

test/curve_analysis/test_guess.py

@@ -186,3 +186,24 @@ def test_baseline_sinusoidal(self, b0, x0, a, freq):
         b0_guess = guess.constant_sinusoidal_offset(y)
 
         self.assertAlmostEqual(b0, b0_guess, delta=0.1)
+
+    @data(
+        # typical 1Q
+        [0.5, 0.5, 0.99],
+        # typical 2Q
+        [0.25, 0.75, 0.97],
+        # alpha around equation switching
+        [0.48, 0.46, 0.85],
+        # bad limit
+        [0.20, 0.36, 0.72],
+        [0.55, 0.40, 0.65],
+    )
+    @unpack
+    def test_rb_decay(self, a, b, alpha):
+        """Test of rb decay basis guess."""
+        x = np.arange(1, 100, 5)
+        y = a * alpha**x + b
+
+        alpha_guess = guess.rb_decay(x, y, a=a, b=b)
+
+        self.assertAlmostEqual(alpha, alpha_guess, delta=alpha * 0.1)