Full, cleaned support for internal Chebyshev basis use, with updated …

…README indicating the change.

README.md

@@ -5,6 +5,10 @@
 
 <!-- *TODO: add some symbols for the main headings, add a title that matches the title of the repo, and overall remove extraneous lists of outputs and inputs outside of a dedicated example section.* -->
 
+## Recent news and updates
+
+> :round_pushpin: The codebase has recently been updated to work internally almost entirely with polynomials in the Chebyshev, rather than the monomial, basis. This improves the numerical stability of the `laurent` method substantially, and should present almost no visible differences for an end-user. Currently the `poly2angles` command-line function is the only major input for monomial basis targets, while the main `QuantumSignalProcessingPhases` method and various helper and plotting methods all accept only Chebyshev basis inputs with internal checks. For conversion between bases, see [below](#computing-qsp-phases-even-faster-with-the-sym-qsp-method).
+
 ## Introduction
 
 [Quantum signal processing](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.118.010501) (QSP) is a flexible quantum algorithm subsuming Hamiltonian simulation, quantum linear systems solvers, amplitude amplification, and many other common quantum algorithms. Moreover, for each of these applications, QSP often exhibits state-of-the-art space and time complexity, with comparatively simple analysis. QSP and its related algorithms, e.g., [Quantum singular value transformation](https://arxiv.org/abs/1806.01838) (QSVT), using basic alternating circuit ansätze originally inspired by composite pulse techniques in NMR, permit one to *transform the spectrum of linear operators by near arbitrary polynomial functions*, with the aforementioned good numerical properties and simple analytic treatment.
@@ -263,9 +267,9 @@ Note that here, unlike in the standard case, we return the full phases (analogou
 > ccoefs = np.polynomial.chebyshev.poly2cheb(pcoefs)
 > pcoefs = np.polynomial.chebyshev.cheb2poly(ccoefs)
 > ```
-> For numerical stability, however, casting between `Polynomial` and `Chebyshev` objects at high degrees is not recommended.
+> For numerical stability, however, casting between `Polynomial` and `Chebyshev` objects at high degrees is not recommended. As such, recent updates have forced almost all methods to work in the Chebyshev basis by default; notably, in our convention, this will also exactly correspond to the Laurent basis up to copying some of the coefficients.
 >
-> As it stands, the `sym_qsp` method works exclusively in the Chebyshev basis, and performs internal checks on passed objects, and numpy support for Chebyshev polynomial optimization and analysis is thorough.
+> As it stands, the `sym_qsp` and most of the `laurent` method now work exclusively in the Chebyshev basis, performing internal checks on passed objects.
 
 Plotting is also supported with the `sym_qsp` method, again calling
 
@@ -348,7 +352,7 @@ usage: pyqsp [options] cmd
 Version: 0.1.6
 Commands:
 
-    poly2angles - compute QSP phase angles for the specified polynomial (use --poly)
+    poly2angles - compute QSP phase angles for the specified polynomial (use --poly). Currently, if using the `laurent` method, this polynomial is expected in monomial basis, while for `sym_qsp` method the Chebyshev basis is expected. Eventually a new flag will be added to specify basis for use with any method.
     hamsim      - compute QSP phase angles for Hamiltonian simulation using the Jacobi-Anger expansion of exp(-i tau sin(2 theta))
     invert      - compute QSP phase angles for matrix inversion, i.e. a polynomial approximation to 1/a, for given delta and epsilon parameter values
     angles      - generate QSP phase angles for the specified --seqname and --seqargs
@@ -376,7 +380,8 @@ Examples:
     pyqsp --plot --func "np.sign(x - 0.5) + np.sign(-1*x - 0.5)" --polydeg 151 --scale 0.9 sym_qsp_func
 
     # Note older examples using the 'laurent' method.
-    pyqsp --plot --tolerance=0.1 --seqargs 2 invert
+    pyqsp --plot --tolerance=0.1 --seqargs 4 invert
+    pyqsp --plot --seqargs=10,0.1 hamsim
     pyqsp --plot-npts=4000 --plot-positive-only --plot-magnitude --plot --seqargs=1000,1.0e-20 --seqname fpsearch angles
     pyqsp --plot-npts=100 --plot-magnitude --plot --seqargs=23 --seqname erf_step angles
     pyqsp --plot-npts=100 --plot-positive-only --plot --seqargs=23 --seqname erf_step angles

pyqsp/angle_sequence.py

@@ -4,17 +4,12 @@
 from numpy.polynomial.polynomial import Polynomial
 from numpy.polynomial.chebyshev import Chebyshev
 
-# TODO: add a Chebyshev import here, which is flip-flopped according to the option given to the main phase-finding method; there should never be a way that we confuse these two paths, which might require entirely deprecating one.
-
 from pyqsp.completion import completion_from_root_finding
 from pyqsp.decomposition import angseq
 from pyqsp.LPoly import LAlg, LPoly
 from pyqsp.poly import StringPolynomial, TargetPolynomial
 from pyqsp.response import ComputeQSPResponse
-
 from pyqsp.sym_qsp_opt import newton_solver
-# TODO: Determine eventually why this import is not wrt pyqsp.
-# import sym_qsp_opt
 
 class AngleFindingError(Exception):
     """Raised when angle finding step fails."""
@@ -75,13 +70,9 @@ def poly2laurent(pcoefs):
     Raises:
         AngleFindingError: if pcoefs does not have definite parity.
     """
-    # convert polynomial coefficients to laurent
-
-    ### TODO: current method to bypass casting to monomial basis.
-
+    # convert polynomial coefficients to laurent; as we are already in the Chebyshev basis, this transformation is trivial
     ccoefs = pcoefs
-
-    ### ccoefs = np.polynomial.chebyshev.poly2cheb(pcoefs)
+    # ccoefs = np.polynomial.chebyshev.poly2cheb(pcoefs) # Conversion in old monomial convention.
 
     # determine parity of polynomial
     is_even = np.max(np.abs(ccoefs[0::2])) > 1e-8
@@ -133,36 +124,34 @@ def QuantumSignalProcessingPhases(
         ValueError: Raised if invalid model (or method) is specified.
     """
 
-    # TODO: we can remove the conditional branch below eventually, as everything should be input as Chebyshev coefficients.
-
-    # Depending on the chebyshev_basis flag, convert coefficient list if provided to Polynomial/Chebyshev object.
-    if not chebyshev_basis:
-        # if isinstance(poly, np.ndarray) or isinstance(poly, list):
-        #     poly = Polynomial(poly)
-        # # For TargetPolynomial class, cast back to parent object.
-        # elif isinstance(poly, Polynomial):
-        #     poly = Polynomial(poly.coef)
-        # else:
-        #     raise ValueError(
-        #         f"Must provide coefficient list/array, or Polynomial object.")
-        ### TODO: current bypass for chebyshev basis
-        if isinstance(poly, np.ndarray) or isinstance(poly, list):
-            poly = Chebyshev(poly)
-        elif isinstance(poly, Chebyshev):
-            pass
-        else:
-            raise ValueError(
-                f"Must provide Chebyshev coefficient list/array, or Chebyshev polynomial object.")
+    # As all inputs are assured to be in the Chebyshev basis, we check whether we have been given a Chebyshev object or an ndarray of coefficients, and cast to Chebyshev object.
+    if isinstance(poly, np.ndarray) or isinstance(poly, list):
+        poly = Chebyshev(poly)
+    elif isinstance(poly, Chebyshev):
+        pass
     else:
-        # TODO: Currently TargetPolynomial is not handled in the Chebyshev basis; we can choose to accommodate this or not.
-        if isinstance(poly, np.ndarray) or isinstance(poly, list):
-            poly = Chebyshev(poly)
-        elif isinstance(poly, Chebyshev):
-            # Currently, there should be no sub-class that we have to cast from in this branch.
-            pass
-        else:
-            raise ValueError(
-                f"Must provide Chebyshev coefficient list/array, or Chebyshev polynomial object.")
+        raise ValueError(
+            f"Must provide Chebyshev coefficient list/array, or Chebyshev polynomial object.")
+
+    # Note: as all polynomials are now guaranteed to be in the Chebyshev basis, we bypass the conditional below.
+    # if not chebyshev_basis:
+    #     if isinstance(poly, np.ndarray) or isinstance(poly, list):
+    #         poly = Polynomial(poly)
+    #     # For TargetPolynomial class, cast back to parent object.
+    #     elif isinstance(poly, Polynomial):
+    #         poly = Polynomial(poly.coef)
+    #     else:
+    #         raise ValueError(
+    #             f"Must provide coefficient list/array, or Polynomial object.")
+    # else:
+    #     if isinstance(poly, np.ndarray) or isinstance(poly, list):
+    #         poly = Chebyshev(poly)
+    #     elif isinstance(poly, Chebyshev):
+    #         # Currently, there should be no sub-class that we have to cast from in this branch.
+    #         pass
+    #     else:
+    #         raise ValueError(
+    #             f"Must provide Chebyshev coefficient list/array, or Chebyshev polynomial object.")
 
     if measurement is None:
         if signal_operator == "Wx":
@@ -193,24 +182,15 @@ def QuantumSignalProcessingPhases(
 
         # Perform completion
         if model in {("Wx", "x"), ("Wz", "z")}:
-            # Capitalization: eps/2 amount of error budget is put to the highest power for sake of numerical stability.
-
-            ### Currently commented for Chebyshev bypass
-            # poly = suc * \
-            #     (poly + Polynomial([0, ] * poly.degree() + [eps / 2, ]))
-
+            # Capitalization: eps/2 amount of error budget is put to the highest power for sake of numerical stability. We are assured to be in the Chebyshev basis.
             poly = suc * \
                  (poly + Chebyshev([0, ] * poly.degree() + [eps / 2, ]))
 
             lcoefs = poly2laurent(poly.coef)
             lalg = completion_from_root_finding(lcoefs, coef_type="F")
         elif model == ("Wx", "z"):
-            # TODO: unlike in the Laurent picture, the 'P' type coefficients are expected in the monomial basis.
-            #
-            # TODO: given the Chebyshev bypass, now is the only setting in which we need to convert back from the Chebyshev basis to the expected monomial basis
-
+            # Unlike in the Laurent picture, the 'P' type coefficients are expected in the monomial basis. Note this is now the only setting in which we need to convert back from the Chebyshev basis to the expected monomial basis, and should be used rarely.
             monomial_coef = np.polynomial.chebyshev.cheb2poly(poly.coef)
-
             lalg = completion_from_root_finding(monomial_coef, coef_type="P")
         else:
             raise ValueError(

pyqsp/completion.py

@@ -126,8 +126,6 @@ def _fg_completion(F, seed):
     """
     # Find all the roots of Id - (p * ~p) that lies within the upper unit
     # circle.
-    #
-    # NOTE: the updated method and the old method both fed a Chebyshev polynomial into this method; when we transform into the Laurent polynomial picture the monomial coefficients are by construction the Chebyshev coefficients of the standard picture.
     poly = (Id - (F * ~F)).coefs
     converted_poly = np.polynomial.Polynomial(poly)
     roots = converted_poly.roots()
@@ -166,6 +164,7 @@ def _fg_completion(F, seed):
     lst = []
     if seed is None:
         seed = np.random.randint(2, size=len(imag_roots) + len(real_roots))
+    # Note: root inversion randomization has been disabled given new stability-ensuring methods implemented elsewhere.
     for i, root in enumerate(imag_roots):
         if 1: # seed[i]: # Commented out as it seems to improve performance.
             root = 1 / root
@@ -183,9 +182,6 @@ def _fg_completion(F, seed):
     gcoefs = np.real(np.fft.fft(coef_fft, 1 << pp))[
         :degree + 1] / (1 << pp)
 
-    # TODO: print(f"NORM: {norm}")
-    # TODO: print(f"GCOEF: {gcoefs[0]}") # Growing too large at high degree.
-
     # Normalization
     G = LPoly(gcoefs * np.sqrt(norm / gcoefs[0]), -len(gcoefs) + 1)
 

pyqsp/main.py

@@ -35,7 +35,7 @@ def CommandLine(args=None, arglist=None):
 Version: {}
 Commands:
 
-    poly2angles - compute QSP phase angles for the specified polynomial (use --poly)
+    poly2angles - compute QSP phase angles for the specified polynomial (use --poly). Currently, if using the `laurent` method, this polynomial is expected in monomial basis, while for `sym_qsp` method the Chebyshev basis is expected. Eventually a new flag will be added to specify basis for use with any method.
     hamsim      - compute QSP phase angles for Hamiltonian simulation using the Jacobi-Anger expansion of exp(-i tau sin(2 theta))
     invert      - compute QSP phase angles for matrix inversion, i.e. a polynomial approximation to 1/a, for given delta and epsilon parameter values
     angles      - generate QSP phase angles for the specified --seqname and --seqargs
@@ -63,7 +63,7 @@ def CommandLine(args=None, arglist=None):
     pyqsp --plot --func "np.sign(x - 0.5) + np.sign(-1*x - 0.5)" --polydeg 151 --scale 0.9 sym_qsp_func
 
     # Note older examples using the 'laurent' method.
-    pyqsp --plot --tolerance=0.1 --seqargs 2 invert
+    pyqsp --plot --tolerance=0.1 --seqargs 4 invert
     pyqsp --plot --seqargs=10,0.1 hamsim
     pyqsp --plot-npts=4000 --plot-positive-only --plot-magnitude --plot --seqargs=1000,1.0e-20 --seqname fpsearch angles
     pyqsp --plot-npts=100 --plot-magnitude --plot --seqargs=23 --seqname erf_step angles
@@ -291,16 +291,16 @@ def float_list(value):
             (phiset, red_phiset, parity) = angle_sequence.QuantumSignalProcessingPhases(
             coefs, **qspp_args)
         else:
-            ### TODO: chebyshev bypass converts everything to Chebyshev; otherwise we should allow the user to specify the basis with a flag
-
+            # QuantumSignalProcessingPhases now expects input in Chebyshev basis in all cases, so we cast (assumed monomial basis) input.
             cheb_coefs = np.polynomial.chebyshev.poly2cheb(coefs)
-
             phiset = angle_sequence.QuantumSignalProcessingPhases(
             cheb_coefs, **qspp_args)
         if args.plot:
+            # Note: while the input to 'poly2angles' is expected in the monomial basis, PlotQSPResponse, like all other internal methods, expects input in the Chebyshev basis.
+            cheb_coefs = np.polynomial.chebyshev.poly2cheb(coefs)
             response.PlotQSPResponse(
                 phiset,
-                pcoefs=coefs,
+                pcoefs=cheb_coefs,
                 signal_operator=args.signal_operator,
                 **plot_args)
 
@@ -362,9 +362,7 @@ def float_list(value):
 
     elif args.cmd == "invert":
         pg = pyqsp.poly.PolyOneOverX()
-
-        # TODO: currently changed so that we always return Chebyshev coefs.
-
+        # Note that in either case, approximating polynomial is returned in Chebyshev basis.
         pcoefs, scale = pg.generate(
             *args.seqargs,
             return_coef=True,
@@ -644,7 +642,7 @@ def float_list(value):
 
         # Note that the polynomial is renormalized once according to its maximum magnitude (both up and down), and then again according to specified scale in a multiplicative way. We ought to be able to specify this, or at least return the metadata.
 
-        # TODO: new addition for manual rescaling within sym_qsp method.
+        # The parameter scale < 1 allows rescaling only within sym_qsp method.
         if args.scale:
             if np.abs(args.scale >= 1) or (args.scale <= 0):
                 raise ValueError(