On Wed, Jan 24, 2024 at 6:29 AM Fabio Matti <somecallmefabio@gmail.com> wrote:
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?us...
This is a tricky sort of decision. The various polynomial functions were designed so that they could be used with different types, Fractions, for instance, and for degrees less than about 50 max. They are sort of a cross between teaching tools and quick prototyping. I agree that for fast production code and big degrees, DCT is the way to go for Chebyshev. There are even some packages out there designed on that basis. I wouldn't complain if someone produced a package with restricted types that was more efficient than what NumPy has, indeed, types are already restricted for curve fitting. NumPy is trending towards specialization these days, I suspect a separate polynomial package would be a better place for such improvements. Chuck