On Wed, Jan 24, 2024 at 6:29 AM Fabio Matti <somecallmefa...@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?usp=sharing
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
_______________________________________________
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
