# # Floating point error propagation in polynomial multiplication with Fast-Fourier Transform # # - *Simple question*: when using the [FFT](https://en.wikipedia.org/wiki/Multiplication_algorithm#Fourier_transform_methods) (or any Fourier transform methods) for multiplying two polynomials, how to deal with floating point error propagation? # # In this [Jupyter notebook](https://www.jupyter.org/), [I](http://perso.crans.org/besson/) will try to understand this phenomena. # ---- # ## Requirements # - The [numpy.polymul](https://docs.scipy.org/doc/numpy/reference/generated/numpy.polymul.html#numpy.polymul) and [numpy.polyfit](https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html#numpy.polyfit) functions. # In: import numpy as np np.version.full_version # ---- # ## Examples # # Let consider the polynomials $P(X) = 1 + X + X^3$ and $Q(X) = X^2 + X^5$, of degrees $n=3$ and $m=5$: # In: P = [ 1, 0, 1, 1] n = len(P) - 1 P, n # In: Q = [1, 0, 0, 1, 0, 0] m = len(Q) - 1 Q, m # Their product is $(PQ)(X)$, of degree $n+m=8$: # \begin{align*}(PQ)(X) &= (1 + X + X^3) (X^2 + X^5) \\ # &= X^2 + X^5 + X^3 + X^7 + X^5 + X^8 \\ # &= X^2 + X^3 + 2 X^5 + X^7 + X^8 # \end{align*} # In: PQ = polymul(P, Q) d = len(PQ) - 1 PQ, d # If we evaluate both $P(X)$ and $Q(X)$, on $n+m$ different points, $\lambda_i \in \mathbb{R}$ or $\in\mathbb{N}$, then we can *fit* a polynomial of degree $n+m = \delta(PQ)$ on these sampling points, and by uniqueness (thanks to [the Fundamental Theorem of Algebra](https://en.wikipedia.org/wiki/Fundamental_theorem_of_algebra#Corollaries)), it will be equal to $(PQ)(X)$. # # The fit can be obtained, for instance, by Lagrange interpolation, which is not so efficient but easy to implement. # Here, I will simply use the [numpy.polyfit](https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html#numpy.polyfit) function. # In: assert d == n + m # ### Naive interpolation values # # Let us consider the points $\lambda_i = i$, $i=0,\dots,n+m$. # In: lambdas = np.arange(0, d + 1) lambdas # In: values_P = np.polyval(P, lambdas) values_P # In: values_Q = np.polyval(Q, lambdas) values_Q # In: values_PQ = values_P * values_Q values_PQ # In: PQ_sampled = np.polyfit(lambdas, values_PQ, d) PQ_sampled # In: PQ # In: np.asarray(np.round(PQ_sampled), dtype=int) # Ok, at least it seems to work! # # But we saw that even with very small degrees ($n=3, m=5$), floating point errors were not so small on these wrongly chosen points $\lambda_i = i$, $i=0,\dots,n+m$. # The largest "should be 0" value (i.e., $\simeq 0$) value was: # In: np.max(np.abs(PQ_sampled)[np.abs(PQ_sampled) < 0.9]) # ### Chebyshev nodes as interpolation values # # Let us consider the points $\lambda_k = \cos\left(\frac{2k-1}{2d} \pi\right)$, $k=1,\dots,1+d=n+m+1$. # These are called the [Chebyshev nodes](https://en.wikipedia.org/wiki/Chebyshev_nodes). # In: lambdas = np.cos(np.pi * (2 * np.arange(1, 2 + d) - 1) / (2 * d)) lambdas # In: values_P = np.polyval(P, lambdas) values_P # In: values_Q = np.polyval(Q, lambdas) values_Q # In: values_PQ = values_P * values_Q values_PQ # In: PQ_sampled2 = np.polyfit(lambdas, values_PQ, d) PQ_sampled2 # In: PQ # In: np.asarray(np.round(PQ_sampled2), dtype=int) # Ok, at least it seems to work! # # But we saw that even with very small degrees ($n=3, m=5$), floating point errors were not so small on these points: the largest "should be 0" value (i.e., $\simeq 0$) value was: # In: np.max(np.abs(PQ_sampled2)[np.abs(PQ_sampled2) < 0.9]) # ### Benchmark # # Stupidly, let us check if our naive implementation of $(P, Q) \mapsto PQ$ by evaluation-and-interpolation is more or less efficient than numpy.polyfit: # In: def mypolymul(P, Q): n = len(P) - 1 m = len(Q) - 1 d = n + m lambdas = np.cos(np.pi * (2 * np.arange(1, 2 + d) - 1) / (2 * d)) values_P = np.polyval(P, lambdas) values_Q = np.polyval(Q, lambdas) values_PQ = values_P * values_Q PQ_sampled = np.polyfit(lambdas, values_PQ, d) # return PQ_sampled return np.asarray(np.round(PQ_sampled), dtype=int) # In: np.polymul(P, Q) # In: mypolymul(P, Q) # In: import warnings warnings.simplefilter('ignore', np.RankWarning) get_ipython().run_line_magic('timeit', 'np.polymul(P, Q)') get_ipython().run_line_magic('timeit', 'mypolymul(P, Q)') # Of course, our implementation is slower. # But on small polynomials, not so slower. # # What about larger polynomials? # In: def random_polynomial(d=10, maxcoef=1): return np.random.randint(low=-maxcoef, high=maxcoef+1, size=d+1) # In: P = random_polynomial() Q = random_polynomial() P, Q get_ipython().run_line_magic('timeit', 'np.polymul(P, Q)') np.polymul(P, Q) get_ipython().run_line_magic('timeit', 'mypolymul(P, Q)') mypolymul(P, Q) assert np.all(np.polymul(P, Q) == mypolymul(P, Q)) # On a larger example: # In: d = 100 maxcoef = 1 get_ipython().run_line_magic('timeit', 'np.polymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))') get_ipython().run_line_magic('timeit', 'mypolymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))') P, Q = random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef) assert np.all(np.polymul(P, Q) == mypolymul(P, Q)) # In: d = 10 maxcoef = 3 get_ipython().run_line_magic('timeit', 'np.polymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))') get_ipython().run_line_magic('timeit', 'mypolymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))') P, Q = random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef) assert np.all(np.polymul(P, Q) == mypolymul(P, Q)) # In: d = 10 maxcoef = 50 get_ipython().run_line_magic('timeit', 'np.polymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))') get_ipython().run_line_magic('timeit', 'mypolymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))') P, Q = random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef) assert np.all(np.polymul(P, Q) == mypolymul(P, Q)) # Our method is slower. # And wrong. # # That's sad. # ---- # ## Conclusion # # Implementing naively the multiplication of 1D polynomials with evaluation-and-interpolation does not work. # # - It is slower that the FFT based method (available in any numerical computation environment), e.g., numpy.polymul in Python with NumPy. # - And it does not work. Booum! # ---- # > Thanks for reading! # # > See [this repo on GitHub](https://github.com/Naereen/notebooks/) for more notebooks, or [on nbviewer.jupyter.org](https://nbviewer.jupyter.org/github/Naereen/notebooks/). # # > That's it for this demo! See you, folks!