A lookup table for fast Python math

December 11, 2008

Numerical programming frequently requires the use of look-up tables. A look-up table is a collection of pre-computed values. When given an “x” value, the table returns a pre-computed “y” value. Look-up tables can be used to speed up numerical codes, when it is faster to look up a value in the table than it is to compute the value. They are also used when the data in the table cannot be computed–for example, experimental data or averaged results from an ensemble of Monte Carlo simulations. Another application is to compute a value when a function cannot be solved algebraically. Assume that you have a formula for a function q(h). You need the value of h for a given value of q, but the formula cannot be algebraically solved to get h(q). Instead, choose a range of h values, compute the function q(h), and store each value in a look-up table. Now you can get h(q) for any value stored in the table. The major limitation of a look-up table is that it cannot return valid results for any value of q which is outside the range of those stored in the table. Depending on its implementation, the table may be able to interpolate to return values between known points.

Scipy includes two handy lookup table classes, but they are sort of hidden. Here is an example of how to create a lookup table using interp1d and UnivariateSpline:

[code language=“Python”]from scipy.interpolate import interp1d, UnivariateSpline from scipy.stats import norm from numpy import arange from numpy.random import uniform import matplotlib.pyplot as plt from time import time

deltax = 0.1 xmin = -5.0 xmax = 5.0

# Set up samples x_samples = arange(xmin, xmax, deltax) pdf_samples = norm.pdf(x_samples) fig = plt.plot(x_samples, pdf_samples, ‘ro’, label=‘Sampled’)

# Interpolate on a finer grid pdf_interp = interp1d(x_samples, pdf_samples, kind=‘linear’) x_interp = arange(xmin+deltax, xmax-deltax, 0.001) plt.plot(x_interp, pdf_interp(x_interp), ‘b-’, label=‘Interpolated’)

# Do the same thing with UnivariateSpline u = UnivariateSpline(x_samples, pdf_samples, k=1, s=0.0) plt.plot(x_interp, u(x_interp), ‘k-’, label=‘Spline’)

plt.xlabel(‘x’) plt.ylabel(‘pdf’) plt.legend() plt.show()

# Time tests start_time = time() for i in range(100000): pdf_interp(uniform(xmin+deltax, xmax-deltax)) print(“interp1d run time: {}".format(time() - start_time))

start_time = time() for i in range(100000): u(uniform(xmin+deltax, xmax-deltax)) print(“UnivariateSpline run time: {}".format(time() - start_time))

x_fine = arange(xmin + deltax, xmax - deltax - 0.001, 1e-6)

start_time = time() pdf_interp(x_fine) print(“interp1d run time: {}".format(time() - start_time))

start_time = time() u(x_fine) print(“UnivariateSpline run time: {}".format(time() - start_time))[/code]

First, create a Numpy array to store the x values for the lookup table. I’m using the pdf of the normal distribution as an example of a function to interpolate. The table uses linear interpolation to compute values between the known points. Higher-order interpolations can be used, but I don’t need them in this case. Higher-order interpolations can introduce distortions for certain types of data, so they are best avoided unless you really understand the data being fitted. Below is a plot illustrating the results of a look-up table:

The red dots are the values stored in the lookup table, the blue line is interpolated using interp1d, and the black line is interpolated using UnivariateSpline (the blue and black lines are indistinguishable). The function was approximated with a lookup table that was created with an evenly spaced array of x values. Notice how the points stored in the table are spaced more closely together in the y dimension for values of x further from the origin. Thus, the approximation becomes less accurate as x approaches 0, which is a consequence of using a linear “x” spacing with a nonlinear function. You could work around this by using a non-uniform array of “x” values to create the lookup table. This is another reason to understand the function you are approximating, instead of blindly placing values into the table.

Peformance

Like all operations in NumPy, interpolation runs much faster if you pass in a Numpy array instead of operating on individual values. On my older Mac laptop, UnivariateSpline is about twice as fast as interp1d when operating on individual floating point numbers, but interp2d is about twice as fast when operating on NumPy arrays:

interp1d run time: 3.5505311489105225
UnivariateSpline run time: 1.5044260025024414
interp1d run time: 0.1413729190826416
UnivariateSpline run time: 0.27219295501708984

Alternatives to SciPy

If you want to understand more about how lookup tables work, or can’t/won’t install SciPy, check out this interpolated lookup table class. I wouldn’t use it for serious math, since it doesn’t use Numpy and is probably not very fast.

Comments

Comment #1 by Optimizing Python code for fast math | as through a mirror dimly

Optimizing Python code for fast math | as through a mirror dimly - Jan 5, 2009

[…] a previous post I described how to use interp1d from scipy to implement a fast lookup table. Well, it turns out that is actually a slow lookup table. The function […]

Comment #2 by Lookup tables and spline fitting in Python | as through a mirror dimly

Lookup tables and spline fitting in Python | as through a mirror dimly - Apr 3, 2009

[…] approximation can be considerably faster. In a previous post, I showed how to use the function interp1d from scipy.interpolate as a lookup table. In a later post, I showed that interp1d is actually rather slow, and […]