Add precision parameter to GQSP (#1299)

* add precision param to GQSP

* autogen notebooks

---------

Co-authored-by: Matthew Harrigan <mpharrigan@google.com>

Author: anurudhp

qualtran/bloqs/hamiltonian_simulation/hamiltonian_simulation_by_gqsp.ipynb

@@ -85,7 +85,7 @@
     "#### Parameters\n",
     " - `walk_operator`: qubitization walk operator of $H$ constructed from SELECT and PREPARE oracles.\n",
     " - `t`: time to simulate the Hamiltonian, i.e. $e^{-iHt}$\n",
-    " - `precision`: the precision $\\epsilon$ to approximate $e^{it\\cos\\theta}$ to a polynomial.\n"
+    " - `precision`: the precision $\\epsilon$ of the final block encoded $e^{-iHt}$. Split into two: half to approximate $e^{it\\cos\\theta}$ to a polynomial, and half to synthesize the underlying GQSP rotations.\n"
    ]
   },
   {

qualtran/bloqs/hamiltonian_simulation/hamiltonian_simulation_by_gqsp.py

@@ -85,7 +85,9 @@ class HamiltonianSimulationByGQSP(GateWithRegisters):
     Args:
         walk_operator: qubitization walk operator of $H$ constructed from SELECT and PREPARE oracles.
         t: time to simulate the Hamiltonian, i.e. $e^{-iHt}$
-        precision: the precision $\epsilon$ to approximate $e^{it\cos\theta}$ to a polynomial.
+        precision: the precision $\epsilon$ of the final block encoded $e^{-iHt}$. Split into two:
+                   half to approximate $e^{it\cos\theta}$ to a polynomial, and half to synthesize
+                   the underlying GQSP rotations.
     """
 
     walk_operator: QubitizationWalkOperator
@@ -108,7 +110,7 @@ def alpha(self):
     @cached_property
     def degree(self) -> SymbolicInt:
         r"""degree of the polynomial to approximate the function e^{it\cos(\theta)}"""
-        return degree_jacobi_anger_approximation(self.t * self.alpha, precision=self.precision)
+        return degree_jacobi_anger_approximation(self.t * self.alpha, precision=self.precision / 2)
 
     @cached_property
     def approx_cos(self) -> Union[NDArray[np.complex128], Shaped]:
@@ -128,6 +130,7 @@ def gqsp(self) -> GeneralizedQSP:
             self.walk_operator,
             self.approx_cos,
             negative_power=self.degree,
+            precision=self.precision / 2,
             verify=True,
             verify_precision=1e-4,
         )

qualtran/bloqs/hamiltonian_simulation/hamiltonian_simulation_by_gqsp_test.py

@@ -107,5 +107,5 @@ def test_hamiltonian_simulation_by_gqsp_t_complexity():
 
     symbolic_hamsim_by_gqsp = _symbolic_hamsim_by_gqsp()
     tau, t, inv_eps = sympy.symbols(r"\tau t \epsilon^{-1}", positive=True)
-    T = big_O(tau * t + sympy.log(inv_eps) / sympy.log(sympy.log(inv_eps)))
+    T = big_O(tau * t + sympy.log(2 * inv_eps) / sympy.log(sympy.log(2 * inv_eps)))
     assert symbolic_hamsim_by_gqsp.t_complexity() == TComplexity(t=T, clifford=T, rotations=T)  # type: ignore[arg-type]

qualtran/bloqs/qsp/generalized_qsp.ipynb

@@ -84,7 +84,8 @@
     " - `U`: Unitary operation.\n",
     " - `P`: Co-efficients of a complex QSP polynomial.\n",
     " - `Q`: Co-efficients of a complex QSP polynomial.\n",
-    " - `negative_power`: value of $k$, which effectively applies $z^{-k} P(z)$. defaults to 0. \n",
+    " - `negative_power`: value of $k$, which effectively applies $z^{-k} P(z)$. defaults to 0.\n",
+    " - `precision`: The error in the synthesized unitary. This is used to compute the required precision for each single qubit SU2 rotation. \n",
     "\n",
     "#### References\n",
     " - [Generalized Quantum Signal Processing](https://arxiv.org/abs/2308.01501). Motlagh and Wiebe. (2023). Theorem 3; Figure 2; Theorem 6.\n"

qualtran/bloqs/qsp/generalized_qsp.py

@@ -31,7 +31,7 @@
 )
 from qualtran.bloqs.basic_gates.su2_rotation import SU2RotationGate
 from qualtran.linalg.polynomial.qsp_testing import assert_is_qsp_polynomial
-from qualtran.symbolics import is_symbolic, Shaped, slen, smax, smin, SymbolicInt
+from qualtran.symbolics import is_symbolic, Shaped, slen, smax, smin, SymbolicFloat, SymbolicInt
 
 if TYPE_CHECKING:
     import cirq
@@ -263,6 +263,8 @@ class GeneralizedQSP(GateWithRegisters):
         P: Co-efficients of a complex QSP polynomial.
         Q: Co-efficients of a complex QSP polynomial.
         negative_power: value of $k$, which effectively applies $z^{-k} P(z)$. defaults to 0.
+        precision: The error in the synthesized unitary. This is used to compute the required
+                   precision for each single qubit SU2 rotation.
 
     References:
         [Generalized Quantum Signal Processing](https://arxiv.org/abs/2308.01501)
@@ -273,6 +275,7 @@ class GeneralizedQSP(GateWithRegisters):
     P: Union[Tuple[complex, ...], Shaped] = field(converter=_to_tuple)
     Q: Union[Tuple[complex, ...], Shaped] = field(converter=_to_tuple)
     negative_power: SymbolicInt = field(default=0, kw_only=True)
+    precision: SymbolicFloat = field(default=1e-11, kw_only=True)
 
     @P.validator
     @Q.validator
@@ -292,6 +295,7 @@ def from_qsp_polynomial(
         P: Union[NDArray[np.number], Sequence[complex], Shaped],
         *,
         negative_power: SymbolicInt = 0,
+        precision: SymbolicFloat = 0,
         verify: bool = False,
         verify_precision=1e-7,
     ) -> 'GeneralizedQSP':
@@ -301,7 +305,7 @@ def from_qsp_polynomial(
         if verify:
             assert_is_qsp_polynomial(P)
         Q = qsp_complementary_polynomial(P, verify=verify, verify_precision=verify_precision)
-        return GeneralizedQSP(U, P, Q, negative_power=negative_power)
+        return GeneralizedQSP(U, P, Q, negative_power=negative_power, precision=precision)
 
     @cached_property
     def _qsp_phases(self) -> Tuple[NDArray[np.floating], NDArray[np.floating], float]:
@@ -311,13 +315,18 @@ def _qsp_phases(self) -> Tuple[NDArray[np.floating], NDArray[np.floating], float
             )
         return qsp_phase_factors(self.P, self.Q)
 
+    @cached_property
+    def _eps_per_rotation(self):
+        """precision to synthesize each SU2 rotation."""
+        return self.precision / (slen(self.P) + 1)
+
     @cached_property
     def signal_rotations(self) -> NDArray[SU2RotationGate]:  # type: ignore[type-var]
         thetas, phis, lambd = self._qsp_phases
 
         return np.array(
             [
-                SU2RotationGate(theta, phi, lambd if i == 0 else 0)
+                SU2RotationGate(theta, phi, lambd if i == 0 else 0, eps=self._eps_per_rotation)
                 for i, (theta, phi) in enumerate(zip(thetas, phis))
             ]
         )