diff --git a/README.md b/README.md index 42478579..3feddbe4 100644 --- a/README.md +++ b/README.md @@ -5,6 +5,10 @@ +## 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 diff --git a/pyqsp/angle_sequence.py b/pyqsp/angle_sequence.py index 1aa2e116..cccdc10a 100644 --- a/pyqsp/angle_sequence.py +++ b/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( diff --git a/pyqsp/completion.py b/pyqsp/completion.py index 46ac842d..65a24558 100644 --- a/pyqsp/completion.py +++ b/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) diff --git a/pyqsp/main.py b/pyqsp/main.py index 0ffca1c8..ee288973 100644 --- a/pyqsp/main.py +++ b/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( diff --git a/pyqsp/poly.py b/pyqsp/poly.py index 5c0968fd..645d940d 100644 --- a/pyqsp/poly.py +++ b/pyqsp/poly.py @@ -39,7 +39,7 @@ def target(self, arg): # ----------------------------------------------------------------------------- -### NOTE: Use of this is curently deprecated, but should be further excised for proper Chebyshev bypass. +#Note: use of this class is deprecated, but retained here for legacy applications. class TargetPolynomial(np.polynomial.Polynomial): ''' Polynomial with ideal target @@ -121,12 +121,12 @@ def generate( print(f"[PolyCosineTX] rescaling by {scale}.") if return_coef: - if chebyshev_basis: - pcoefs = g.coef - else: - ### TODO: new bypass for Chebyshev coefficients. - pcoefs = g.coef - # pcoefs = np.polynomial.chebyshev.cheb2poly(g.coef) + pcoefs = g.coef + # Note that now all approximating polynomials are in the Chebyshev basis, bypassing the branching conditional below. + # if chebyshev_basis: + # pcoefs = g.coef + # else: + # pcoefs = np.polynomial.chebyshev.cheb2poly(g.coef) if ensure_bounded and return_scale: return pcoefs, scale else: @@ -177,12 +177,12 @@ def generate( print(f"[PolySineTX] rescaling by {scale}.") if return_coef: - if chebyshev_basis: - pcoefs = g.coef - else: - ### TODO: new bypass for Chebyshev coefficients. - pcoefs = g.coef - # pcoefs = np.polynomial.chebyshev.cheb2poly(g.coef) + pcoefs = g.coef + # Note: all polynomials are expected in the Chebyshev basis, bypassing the conditional below. + # if chebyshev_basis: + # pcoefs = g.coef + # else: + # pcoefs = np.polynomial.chebyshev.cheb2poly(g.coef) if ensure_bounded and return_scale: return pcoefs, scale else: @@ -265,7 +265,7 @@ def generate( print("[PolyOneOverX] bounding to 0.5") g = scale * g - # # Internal plotting for debugging purposes. + # TODO: Internal plotting for debugging purposes. # npts = 500 # adat = np.linspace(-1, 1, npts) # plt.plot(adat, g(adat)) @@ -276,14 +276,12 @@ def generate( # plt.show() if return_coef: - if chebyshev_basis: - pcoefs = g.coef - else: - ### TODO: current change to check what happens for chebyshev - - pcoefs = g.coef - - ### pcoefs = np.polynomial.chebyshev.cheb2poly(g.coef) + pcoefs = g.coef + # Note: as all polynomials are now treated in the Chebyshev basis, the conditional below has been bypassed. + # if chebyshev_basis: + # pcoefs = g.coef + # else: + # pcoefs = np.polynomial.chebyshev.cheb2poly(g.coef) # print(f"[pyqsp.PolyOneOverX] pcoefs={pcoefs}") if ensure_bounded and return_scale: return pcoefs, scale @@ -294,8 +292,6 @@ def generate( # ----------------------------------------------------------------------------- -# TODO: for the moment this is causing PolyRect to return a polynomial with NaN coefficients for half the terms and trivial coefficients for the others. - class PolyOneOverXRect(PolyGenerator): def help(self): @@ -320,25 +316,30 @@ def generate( coefs_rect, scale2 = PolyRect().generate(degree, delta, kappa, - epsilon, # TODO: Added to avoid scaling error. + epsilon, ensure_bounded, return_scale=True, chebyshev_basis=chebyshev_basis) - # Modified to handle Chebyshev basis for convolution. - # - ### TODO: Chebyshev bypass - if 0: # not chebyshev_basis: - poly_invert = np.polynomial.Polynomial(coefs_invert) - poly_rect = np.polynomial.Polynomial(coefs_rect) - pcoefs = (poly_invert * poly_rect).coef - else: - poly_invert = np.polynomial.chebyshev.Chebyshev(coefs_invert) - poly_rect = np.polynomial.chebyshev.Chebyshev(coefs_rect) + poly_invert = np.polynomial.chebyshev.Chebyshev(coefs_invert) + poly_rect = np.polynomial.chebyshev.Chebyshev(coefs_rect) + + mult_result = np.polynomial.chebyshev.chebmul(poly_invert.coef, poly_rect.coef) + pcoefs = mult_result + + # Note: all polynomials are now considered in Chebyshev basis, bypassing the original conditional below. + # if not chebyshev_basis: + # poly_invert = np.polynomial.Polynomial(coefs_invert) + # poly_rect = np.polynomial.Polynomial(coefs_rect) - mult_result = np.polynomial.chebyshev.chebmul(poly_invert.coef, poly_rect.coef) - pcoefs = mult_result + # pcoefs = (poly_invert * poly_rect).coef + # else: + # poly_invert = np.polynomial.chebyshev.Chebyshev(coefs_invert) + # poly_rect = np.polynomial.chebyshev.Chebyshev(coefs_rect) + + # mult_result = np.polynomial.chebyshev.chebmul(poly_invert.coef, poly_rect.coef) + # pcoefs = mult_result if return_scale: return pcoefs, scale1 * scale2 @@ -381,8 +382,7 @@ def taylor_series( We also evaluate the mean absolute difference on equispaced points over the interval [-1,1]. ''' - ### TOOD: Current chebyshev bypass. - + # Note: PolyTaylorSeries now no longer generates approximating Taylor polynomials, but only Chebyshev interpolations as contained in the assured branch indicated below. This exhibits better stability and convergence. if 1: # chebyshev_basis: cheb_samples = 2*degree # Set to prevent aliasing; note that all methods calling TaylorSeries implicitly have their cheb_samples specifications overruled here. # Generate x and y values for fit according to func; note use of chebyshev nodes of the first kind. diff --git a/pyqsp/response.py b/pyqsp/response.py index 3a61df72..1488b1f0 100644 --- a/pyqsp/response.py +++ b/pyqsp/response.py @@ -167,10 +167,7 @@ def PlotQSPResponse( if pcoefs is not None: poly = pcoefs - ### TODO: Chebyhsve plotting bypass - # if not isinstance(poly, np.polynomial.Polynomial): - # poly = np.polynomial.Polynomial(pcoefs) - # NOTE: we cast instead to the expected Chebyshev type. + # Note: all polynomials are expected in the Chebyshev basis, and so we case the target polynomial accordingly. if not isinstance(poly, np.polynomial.chebyshev.Chebyshev): poly = np.polynomial.chebyshev.Chebyshev(pcoefs) expected = poly(adat) diff --git a/pyqsp/sym_qsp_opt.py b/pyqsp/sym_qsp_opt.py index 6d131440..e78278a7 100644 --- a/pyqsp/sym_qsp_opt.py +++ b/pyqsp/sym_qsp_opt.py @@ -240,9 +240,7 @@ def newton_solver(coef, parity, **kwargs): if 'crit' in kwargs: crit = kwargs['crit'] - # crit = 5e-2 # TODO: remove after testing. else: - # crit = 5e-2 # TODO: replace with crit = 1e-12. crit = 1e-12 if 'maxiter' in kwargs: maxiter = kwargs['maxiter'] diff --git a/pyqsp/sym_qsp_plotting.py b/pyqsp/sym_qsp_plotting.py index e92c3cc5..6ae4853b 100644 --- a/pyqsp/sym_qsp_plotting.py +++ b/pyqsp/sym_qsp_plotting.py @@ -252,7 +252,7 @@ def main(): freq = 16 pg = poly.PolySign() - # TODO: note that definition of PolySign has been changed to return bare pcoefs and not TargetPolynomial + # Note that definition of PolySign has been changed to return bare pcoefs and not TargetPolynomial, which is now deprecated. pcoefs, scale = pg.generate( degree=161, delta=25, diff --git a/pyqsp/test/test_poly.py b/pyqsp/test/test_poly.py index 4dbfe1ad..2c59af37 100644 --- a/pyqsp/test/test_poly.py +++ b/pyqsp/test/test_poly.py @@ -40,7 +40,7 @@ def test_oneoverx1(self): assert diff < 0.1 # Currently silenced given instability of high-degree completion. - def _test_poly_one_over_x_response1(self): + def test_poly_one_over_x_response1(self): pg = pyqsp.poly.PolyOneOverX() pcoefs = pg.generate(3, 0.3, return_coef=True, ensure_bounded=True) phiset = angle_sequence.QuantumSignalProcessingPhases(