Hey Everybody,
I am approximating the derivative of nonperiodic functions on [-1,1]
using Chebyshev polynomials. The implementation is straightforward and
works well but is painfully slow. The routine wastes most of its time on
a trivial operation (see comment in the code)
Unfortunately the spectral coefficients are calculated recursively and
thus I haven't found a way to speed up the code. I am aware of list
comprehensions, the speed advantage of calling native numpy functions
etc. But no luck yet finding an alternate formulation that speeds up the
calculation. I attached a stripped down variant of the code. I am very
grateful for tips or directions where to look for a solution (well I
know writing a small C extension might be the way to go but Id love to
have speedy Python)
Note: The DCT takes almost no time compared to the statement in the loop.
Cheers
Robert
import scipy.fftpack
import numpy as np
import timeit
def chebyshev_derivative_gc():
""" Calculate the first n derivatives of the function given at the discrete
data points y.
:param y:
:param n:
"""
y = np.random.random( ( 6, 64 ) )
N = 64
axis = -1
n = 2
a = scipy.fftpack.dct( y, type = 2, norm = None, axis = axis )
# Storage array to hold the result
grad = np.zeros( ( n, ) + y.shape )
# TODO: This loop should go into some type of C extension
for i in xrange( n ):
b = np.zeros( y.shape )
b[:, N - 2] = 2.0 * ( N - 1 ) * a[:, N - 1]
# Now calculate coefficients for the first derivative
# TODO: This is the hot loop that wastes 90% of the computing time
for j in xrange( N - 2, 0, -1 ):
# The next line is the real offender - try by commenting out.
b[:, j - 1] = ( 2.0 * j ) * a[:, j] + b[:, j + 1]
pass
a = b.copy()
grad[i] = scipy.fftpack.dct( b, type = 3, norm = None, axis = axis ) / ( 2 * N )
return grad.real
# Timing code
t = timeit.Timer( chebyshev_derivative_gc )
print t.repeat( 3, 500 )
_______________________________________________
NumPy-Discussion mailing list
[email protected]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
