In this Jupyter notebook, I will try to understand this phenomena.
import numpy as np
np.version.full_version
Let consider the polynomials $P(X) = 1 + X + X^3$ and $Q(X) = X^2 + X^5$, of degrees $n=3$ and $m=5$:
P = [ 1, 0, 1, 1]
n = len(P) - 1
P, n
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*} $$
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), 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 function.
assert d == n + m
Let us consider the points $\lambda_i = i$, $i=0,\dots,n+m$.
lambdas = np.arange(0, d + 1)
lambdas
values_P = np.polyval(P, lambdas)
values_P
values_Q = np.polyval(Q, lambdas)
values_Q
values_PQ = values_P * values_Q
values_PQ
PQ_sampled = np.polyfit(lambdas, values_PQ, d)
PQ_sampled
PQ
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:
np.max(np.abs(PQ_sampled)[np.abs(PQ_sampled) < 0.9])
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.
lambdas = np.cos(np.pi * (2 * np.arange(1, 2 + d) - 1) / (2 * d))
lambdas
values_P = np.polyval(P, lambdas)
values_P
values_Q = np.polyval(Q, lambdas)
values_Q
values_PQ = values_P * values_Q
values_PQ
PQ_sampled2 = np.polyfit(lambdas, values_PQ, d)
PQ_sampled2
PQ
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:
np.max(np.abs(PQ_sampled2)[np.abs(PQ_sampled2) < 0.9])
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:
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)
np.polymul(P, Q)
mypolymul(P, Q)
import warnings
warnings.simplefilter('ignore', np.RankWarning)
%timeit np.polymul(P, Q)
%timeit mypolymul(P, Q)
Of course, our implementation is slower. But on small polynomials, not so slower.
What about larger polynomials?
def random_polynomial(d=10, maxcoef=1):
return np.random.randint(low=-maxcoef, high=maxcoef+1, size=d+1)
P = random_polynomial()
Q = random_polynomial()
P, Q
%timeit np.polymul(P, Q)
np.polymul(P, Q)
%timeit mypolymul(P, Q)
mypolymul(P, Q)
assert np.all(np.polymul(P, Q) == mypolymul(P, Q))
On a larger example:
d = 100
maxcoef = 1
%timeit np.polymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))
%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))
d = 10
maxcoef = 3
%timeit np.polymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))
%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))
d = 10
maxcoef = 50
%timeit np.polymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))
%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.
Implementing naively the multiplication of 1D polynomials with evaluation-and-interpolation does not work.
numpy.polymul in Python with NumPy.Thanks for reading!
See this repo on GitHub for more notebooks, or on nbviewer.jupyter.org.
That's it for this demo! See you, folks!