Adds new Bayesian Inference based factory, BayesFactory (unitaryfund#471

)

* Adds BayesFactory

Adds a new factory based that uses Bayesian inference to extrapolate to the zero-noise-limit.

* Adds tests and fixes docstring

* Lint

* Fix BayesFactory test

* Update BayesFactory and CHANGELOG.md

* Change target_accept

* Lint

* Rename to ExpBayesFactory and update docstrings

* Minor fixes

* Remove docstring in reduce

* Clarify noise parameter and fix conflicts

* Update Bayesian Model

Bayesian estimate of the zero-noise limit

* Update requirements.txt

* format style with black

* Update requirements.txt

* Update test_inference.py

* Update inference.py

* Lint

* Attemping to fix some of the CI errors

* flake and black should work now?

* Fixing sphinx error with blank lines

Co-authored-by: Nathan Shammah <[email protected]>
Co-authored-by: Sarah Kaiser <[email protected]>

CHANGELOG.md

@@ -63,6 +63,7 @@ in the *Mitiq Examples* section.
 - Add a new FakeNodesFactory class based on an alternative interpolation method (@elmandouh, gh-444).
 - Add parameter calibration method to find base noise for parameter noise (@yhindy, gh-411).
 - Remove duplication of the reduce method in every (non-adaptive) factory (@elmandouh, gh-470).
+- Add a new ExpBayesFactory class based on Bayesian Inference (@elmandouh, gh-471).
 
 ## Version 0.4.1 (January 12th, 2021)
 

mitiq/zne/inference.py

@@ -34,6 +34,8 @@
 from numpy.lib.polynomial import RankWarning
 from scipy.optimize import curve_fit, OptimizeWarning
 
+import pymc3 as pm
+
 from mitiq import QPROGRAM
 from mitiq.collector import Collector
 
@@ -1479,6 +1481,146 @@ def _zne_curve(scale_factor: float) -> float:
 ]
 
 
+class ExpBayesFactory(BatchedFactory):
+    """Factory object implementing a zero-noise extrapolation algorithm based on
+    Bayesian Inference. The exponential ansatz for the expectation value:
+    E(lambda) = a + b * e**(-c*lambda),
+    where a, b and c are model parameters that need to be estimated,
+    while lambda is the noise scale factor.
+
+    Args:
+        scale_factors: Sequence of noise scale factors at which
+            expectation values should be measured.
+        shot_list: Optional sequence of integers corresponding to the number
+            of samples taken for each expectation value. If this
+            argument is explicitly passed to the factory, it must have
+            the same length of scale_factors and the executor function
+            must accept "shots" as a valid keyword argument.
+
+    Raises:
+        ValueError: If data is not consistent with the extrapolation model.
+    Note:
+        RichardsonFactory and LinearFactory are special cases of PolyFactory.
+    """
+
+    @staticmethod
+    def _exp_ansatz(
+        a: float, b: float, c: float, scale_factor: float
+    ) -> float:
+        """
+        Calculates the expectation given a scale factor and model
+        parameters.
+
+        Args:
+            a: Model parameter.
+            b: Model parameter.
+            c: Model parameter.
+            scale_factors: The array of noise scale factors.
+
+        Returns:
+            The expected value according to the model parameters and
+            scale factor.
+        """
+
+        return a + b * np.exp(-1.0 * c * scale_factor)
+
+    @staticmethod
+    def extrapolate(
+        scale_factors: Sequence[float],
+        exp_values: Sequence[float],
+        full_output: bool = False,
+    ) -> Union[
+        float,
+    ]:
+        """Static method which evaluates an exponential extrapolation to the
+        zero-noise limit using Bayesian inference.
+
+        Args:
+            scale_factors: The array of noise scale factors.
+            exp_values: The array of expectation values.
+            full_output: If False (default), only the zero-noise limit is
+                returned. If True, additional information about the
+                extrapolated limit is returned too.
+
+        Returns:
+            zne_limit: The extrapolated zero-noise limit. If "full_output" is\
+                False (default value), only this parameter is returned.
+            opt_params: The parameter array of the best fitting model.
+            zne_curve: The callable function which best fit the input data.\
+                It maps a real noise scale factor to a real expectation value.\
+                It is equal "zne_limit" when evaluated at zero.
+
+        Note:
+            This method computes the zero-noise limit only from the information
+            contained in the input arguments. To extrapolate from the internal
+            data of an instantiated Factory object, the bound method
+            ".reduce()" should be called instead.
+        """
+
+        with pm.Model():
+            """
+            We assume that the priors for the model parameters are normally
+            distributed, while the standard deviation is half normal.
+            """
+            a = pm.Normal("a", 0, 5)
+            b = pm.Normal("b", 0, 5)
+            c = pm.Normal("c", 0, 5)
+            sd = pm.HalfNormal("sd", 10)
+
+            pm.Normal(
+                "expval",
+                mu=ExpBayesFactory._exp_ansatz(a, b, c, scale_factors),
+                sd=sd,
+                observed=exp_values,
+            )
+
+            trace = pm.sample(target_accept=0.95, return_inferencedata=False)
+
+        def zne_curve(scale_factor: float) -> float:
+            samples = []
+            n_samples = len(trace["a"])
+
+            for i in range(1000, n_samples):
+                sample = ExpBayesFactory._exp_ansatz(
+                    trace["a"][i], trace["b"][i], trace["c"][i], scale_factor,
+                )
+
+                samples.append(sample)
+
+            return np.mean(samples)
+
+        # Optimal model parameters:
+        a = trace["a"].mean()
+        b = trace["b"].mean()
+        c = trace["c"].mean()
+        opt_params, zne_error = [a, b, c], trace["sd"].mean()
+
+        zne_limit = zne_curve(0)
+        if not full_output:
+            return zne_limit
+
+        params_cov = None
+        return zne_limit, zne_error, opt_params, params_cov, zne_curve
+
+    def reduce(self) -> float:
+        (
+            self._zne_limit,
+            self._zne_error,
+            self._opt_params,
+            self._params_cov,
+            self._zne_curve,
+        ) = self.extrapolate(  # type: ignore
+            self.get_scale_factors(),
+            self.get_expectation_values(),
+            full_output=True,
+        )
+        self._already_reduced = True
+        return self._zne_limit
+
+    def __eq__(self, other: Any) -> bool:
+        return BatchedFactory.__eq__(self, other)
+
+
 class AdaExpFactory(AdaptiveFactory):
     """Factory object implementing an adaptive zero-noise extrapolation
     algorithm assuming an exponential ansatz y(x) = a + b * exp(-c * x),

mitiq/zne/tests/test_inference.py

@@ -34,6 +34,7 @@
     ExpFactory,
     PolyExpFactory,
     AdaExpFactory,
+    ExpBayesFactory,
 )
 
 
@@ -126,6 +127,13 @@ def f_runge(x: float) -> float:
     return 1.0 / ((x - 2) ** 2 + 1.0)
 
 
+def f_exp_bayes(
+    x: float, err: float = STAT_NOISE, rnd_state: RandomState = np.random
+) -> float:
+    """Exponential decay with with parameters suiting BayesFactroy."""
+    return 0.5 + 0.5 * np.exp(-0.5 * x) + rnd_state.normal(scale=err)
+
+
 @mark.parametrize("test_f", [f_lin, f_non_lin])
 def test_noise_seeding(test_f: Callable[[float], float]):
     """Check that seeding works as expected."""
@@ -147,6 +155,7 @@ def test_noise_seeding(test_f: Callable[[float], float]):
         PolyFactory,
         ExpFactory,
         PolyExpFactory,
+        ExpBayesFactory,
     ),
 )
 def test_get_scale_factors_static_factories(factory):
@@ -204,6 +213,7 @@ def test_get_scale_factors_adaptive_factories(factory):
         PolyFactory,
         ExpFactory,
         PolyExpFactory,
+        ExpBayesFactory,
     ),
 )
 def test_get_expectation_values_static_factories(factory):
@@ -273,6 +283,7 @@ def test_get_expectation_values_adaptive_factories(factory):
         PolyFactory,
         ExpFactory,
         PolyExpFactory,
+        ExpBayesFactory,
     ),
 )
 @mark.parametrize("batched", (True, False))
@@ -316,6 +327,7 @@ def executor(circuit):
         PolyFactory,
         ExpFactory,
         PolyExpFactory,
+        ExpBayesFactory,
     ),
 )
 def test_run_batched_with_keyword_args_list(factory):
@@ -440,6 +452,18 @@ def test_poly_extr():
     )
 
 
+def test_exp_bayes_extr():
+    """Test of the ExpBayesFactory's extrapolator."""
+    x_vals = np.linspace(1.0, 5.0, 20)
+    seeded_f = apply_seed_to_func(f_exp_bayes, SEED)
+    fac = ExpBayesFactory(scale_factors=x_vals)
+    assert not fac._opt_params
+    fac.run_classical(seeded_f)
+    zne_value = fac.reduce()
+    assert np.isclose(zne_value, seeded_f(0, err=0), atol=CLOSE_TOL)
+    assert np.isclose(fac._zne_curve(0), seeded_f(0, err=0), atol=CLOSE_TOL)
+
+
 @mark.parametrize("order", [2, 3, 4, 5])
 def test_opt_params_poly_factory(order):
     """Tests that optimal parameters are stored after calling the reduce
@@ -646,7 +670,8 @@ def test_avoid_log_keyword():
 
 
 @mark.parametrize(
-    "factory", (LinearFactory, RichardsonFactory, FakeNodesFactory)
+    "factory",
+    (LinearFactory, RichardsonFactory, FakeNodesFactory, ExpBayesFactory,),
 )
 def test_too_few_scale_factors(factory):
     """Test less than 2 scale_factors."""
@@ -850,7 +875,8 @@ def test_params_cov_and_zne_std():
 
 
 @mark.parametrize(
-    "factory", [LinearFactory, RichardsonFactory, FakeNodesFactory]
+    "factory",
+    [LinearFactory, RichardsonFactory, FakeNodesFactory, ExpBayesFactory,],
 )
 def test_execute_with_zne_fit_fail(factory):
     """Tests errors are raised when asking for fitting parameters that can't

requirements.txt

@@ -1,3 +1,4 @@
 numpy~=1.20.1
 scipy~=1.4.1
 cirq~=0.10.0
+pymc3~=3.11.2