Hi,
In the `numpy.polynomial.chebyshev` module, the function for raising a
Chebyshev polynomial to a power, `chebpow` [1], is essentially implemented in
the following way:
{{{#!highlight python
def chebpow(c, pow):
"""Raise a Chebyshev series to a power."""
zs = _cseries_to_zseries(c)
prd = zs
for i in range(2, pow + 1):
prd = np.convolve(prd, zs)
return _zseries_to_cseries(prd)
}}}
For large coefficient arrays `c` and big exponents `pow`, this procedure is not
efficient. In fact, the complexity of this function is `O(pow*len(c)^2)`, since
the numpy convolution does not make use of a Fast Fourier Transform (FFT).
It is known that Chebyshev polynomials can be multiplied with Discrete Cosine
Transforms (DCT) [2]. What results is the following
algorithm`O(pow*len(c)*log(pow*len(c)))` algorithm for raising a Chebyshev
polynomial with coefficients `c` to an integer power:
{{{#!highlight python
def chebpow_dct(c, pow):
"""Raise a Chebyshev series to a power."""
pad_length = (pow - 1) * (len(c) - 1)
c = np.pad(c, (0, pad_length))
c[1:-1] /= 2
c_pow = idct(dct(c) ** pow)
c_pow[1:-1] *= 2
return c_pow
}}}
The only issue I am having is that as far as I know, `numpy` (unlike `scipy`)
does not have a specialized implementation for the DCT. So the only way of
getting the code to work is "emulating" a DCT with two calls to
`numpy.fft.rfft`, which is slightly slower than using the `scipy.fft.dct`.
I have created a Google colab notebook which compares the error and runtime of
the different implementations (current implementation, implementation using
`scipy.fft.dct`, and pure numpy implementation) [3]. Especially for larger
degree polynomials and higher powers this enhancement would make a huge
difference in terms of runtime.
Similarly, `chebmul` and `chebinterpolate` can also be implemented more
efficiently by using a DCT.
Do you think this enhancement is worth pursuing, and should I create a
pull-request for it?
Best,
Fabio
[1] https://github.com/numpy/numpy/blob/main/numpy/polynomial/chebyshev.py#L817
[2] https://www.sciencedirect.com/science/article/pii/0024379595006966
[3]
https://colab.research.google.com/drive/1JtDDeWC1CEQHDidZ9f5_Ma_ifoBv4Tuz?usp=sharing
_______________________________________________
NumPy-Discussion mailing list -- numpy-discussion@python.org
To unsubscribe send an email to numpy-discussion-le...@python.org
https://mail.python.org/mailman3/lists/numpy-discussion.python.org/
Member address: arch...@mail-archive.com
