Computing complementary polynomials¶
An important step in the synthesis of a QSP protocol for a desired polynomial transformation \(P(z)\) is the completion problem: compute a complementary polynomial, i.e., a polynomial \(Q(z)\) such that \(|P(z)|^2 + |Q(z)|^2 = 1\) whenever \(|z| = 1\).
Using Weiss' algorithm recommended¶
nlft-qsp provides the weiss module which takes care of this.
from nlft_qsp import *
from nlft_qsp.solvers import weiss
P = Polynomial([1/3, 1/5, 1/7])
Q = weiss.complete(P)
print("P(z) =", P)
print("Q(z) =", Q)
print("|P|^2 + |Q|^2 = ", P * P.conjugate() + Q * Q.conjugate())
P(z) = (0.3333333333333333+0j) z^0 + (0.2+0j) z^1 + (0.14285714285714285+0j) z^2
Q(z) = (-0.052808906981778855-2.5762378795436724e-18j) z^-2 + (-0.11218803842422881+3.2738629385311154e-18j) z^-1 + (0.9017237875322446+4.595282034157402e-19j) z^0
|P|^2 + |Q|^2 = (-1.3877787807814457e-17-3.851859888774472e-35j) z^-2 + (-1.1102230246251565e-16-3.851859888774472e-35j) z^-1 + (1.0000000000000002-3.851859888774472e-35j) z^0 + (5.551115123125783e-17-3.851859888774472e-35j) z^1 + (4.85722573273506e-17-3.851859888774472e-35j) z^2
Note
The \(Q\) computed above actually has non-positive frequencies. This because weiss.complete() actually returns a complementary polynomial following the convention for the nonlinear Fourier transform (see the definition here).
It is possible to obtain a normal polynomial with non-negative frequencies by using Polynomial.shift().
Computing the ratio \(P/Q\)¶
First approaches based on nonlinear Fourier analysis (namely riemann_hilbert and half_cholesky) do not need the complementary polynomial \(Q\), but the Laurent expansion of the ratio \(P/Q\). For this reason, the weiss module also conveniently provides a ratio() method for this. Note that it also returns \(Q\) itself.
from nlft_qsp import *
from nlft_qsp.solvers import weiss
P = Polynomial([1/3, 1/5, 1/7])
Q, R = weiss.ratio(P)
print("P(1)/Q(1) =", P(1)/Q(1))
print("R(1) =", R(1))
Using Prony's method¶
The package also provides the prony module (taken and adapted from this repo), which completes polynomials based on Prony's method.
from nlft_qsp import *
from nlft_qsp.solvers import prony
P = Polynomial([1/3, 1/5, 1/7])
Q = prony.complete(P)
print("P(z) =", P)
print("Q(z) =", Q)
print("|P|^2 + |Q|^2 = ", P * P.conjugate() + Q * Q.conjugate())
P(z) = (0.3333333333333333+0j) z^0 + (0.2+0j) z^1 + (0.14285714285714285+0j) z^2
Q(z) = (-0.05280890698177941+3.201254854395866e-17j) z^-2 + (-0.11218803842422782+5.8292190703202056e-18j) z^-1 + (0.9017237875322445+2.1246486134915392e-17j) z^0
|P|^2 + |Q|^2 = (-5.273559366969494e-16+0j) z^-2 + (8.326672684688674e-16+6.5257272063021025e-18j) z^-1 + (0.9999999999999999-1.0558848421106699e-17j) z^0 + (9.43689570931383e-16+1.0558848421106699e-17j) z^1 + (-5.273559366969494e-16-6.5257272063021025e-18j) z^2
Warning
This is provided only for comparison purposes and should not be used, as it does not guarantee numerical stability for large polynomial degrees.
From the Command Line Interface¶
It is possible to compute the complementary polynomial using qspx complete.
The support of \(Q\) will follow the support of \(P\), i.e., if \(P\) has frequencies in \([a, b]\) then also \(Q\) will be shifted to have frequencies in the same range. This is to guarantee that \((P, Q)\) is the actual pair of polynomials to be implementable by QSP.
By passing --nlft the NLFT convention will be chosen, i.e., the support of \(Q\) will end at the zero frequency, in such a way that the pair \((Q, P)\) is in the image of the nonlinear Fourier transform. See here for more information about the nonlinear Fourier transform over \(SU(2)\).