Numerical Python
For the past few months, I've been covering different software packages for scientific computations. For my next several articles, I'm going to be focusing on using Python to come up with your own algorithms for your scientific problems. Python seems to be completely taking over the scientific communities for developing new code, so it is a good idea to have a solid working knowledge of the fundamentals so you can build solutions to your own problems.
In this article, I start with NumPY. My next article will cover SciPy, and then I'll look at some of the more complicated modules available in the following article.
So, let's start with the root module from which almost all other scientific modules are built, NumPY. Out of the box, Python supports real numbers and integers. You also can create complicated data structures with lists, sets and so on. This makes it very easy to write algorithms to solve scientific problems. But, just diving in naively, without paying attention to what is happening under the hood, will lead to inefficient code. This is true with all programming languages, not just Python. Most scientific code needs to squeeze every last available cycle out of the hardware. One of the things to remember about Python is that it is a dynamic language where almost all functions and operators are polymorphic. This means that Python doesn't really know what needs to be done, at a hardware level, until it hits that operation. Unfortunately, this rules out any optimizations that can be made by rearranging operations to take advantage of how they are stored in memory and cache.
One property of Python that causes a much larger problem is polymorphism. In this case, Python needs to check the operands of any operator or function to see what type it is, decide whether this particular operand or function can handle these data types, then use the correct form of the operand or function to do the actual operation. In most cases, this is not really an issue because modern computers have become so fast. But in many scientific algorithms, you end up applying the same operations to thousands, or millions, of data points. A simple example is just taking the square of the first 100,000 numbers:
import time
a = range(100000)
c = []
starttime = time.clock()
for b in a:
c.append(b*b)
endtime = time.clock()
print "Total time for loop: ", (endtime  starttime)
This little program uses the range
function to create a list of the
first 100,000 integers. You need to import the time
module to get access
to the system clock to do timing measurements. Running this on my desktop
takes approximately 0.037775 seconds (remember always to make several
measurements, and take the average). It takes this much time because
for every iteration of the loop, Python needs to check the type of the
b variable before trying to use the multiplication operator. What can
you do if you have even more complicated algorithms?
This is where NumPY comes in. The key element that NumPY introduces is an Ndimensional array object. The great flexibility of Python lists, allowing all sorts of different types of elements, comes at a computational cost. NumPY arrays deal with this cost by introducing restrictions. Arrays can be multidimensional, and they must all be the same data type. Once this is done, you can start to take some shortcuts by taking advantage of the fact that you already know what the type of the elements is. It adds extra overloading functions for the common operators and functions to help optimize uses of arrays.
All of the
normal arithmetic operators work on NumPY arrays in an elementwise
fashion. So, to get the squares of the elements in an array, it would
be as simple as array1 * array1
.
NumPY also uses external standard, optimized libraries written in C or FORTRAN to handle many of the actual manipulations on these array data types. This is handled by libraries like BLAS or lapack. Python simply does an initial check of each of the arrays in question, then hands them as a single object to the external library. The external library does all of the hard work, then hands back a single object containing the result. This removes the need for Python to check each element when using the loop notation above. Using NumPY, the earlier example becomes:
import numpy
import time
a = numpy.arange(1000000)
starttime = time.clock()
c = a * a
endtime = time.clock()
print "Total time used: ", (endtime  starttime)
Running this toy code results in an average run time of 0.011167 seconds for this squaring operation. So the time is cut by one third, and the readability of the code is simplified by getting rid of the loop construct.
I've dealt only with onedimensional arrays so far, but NumPY supports multidimensional arrays just as easily. If you want to define a twodimensional array, or a matrix, you can set up a small one with something like this:
a = numpy.array([[1,2,3,4], [2,3,4,5]])
Basically, you are creating a list of lists, where each of the sublists
is each of the rows of your matrix. This will work only if each of the
sublists is the same size—that is, each of the rows has the same number
of columns. You can get the total number of elements in this matrix,
with the property a.size
. The dimensions of the matrix are stored in
the property a.shape
. In this case, the size is 8, and the shape is
(2, 4), or two rows and four columns. What shape did the array in the
earlier example have? If you go ahead and check, you should see that
the shape is (1000000). The other key properties of these arrays are:

ndim: the number of dimensions.

dtype: the data type of the elements.

itemsize: the size in bytes of each element.

data: the buffer that stores the actual data.
You also can reshape arrays. So if you wanted to take the earlier example of the first 100,000 integers and turn it into a threedimensional array, you could do something like this:
old_array = numpy.arange(100000)
new_array = old_array.reshape(10,100,100)
This will give you a new 3D array laid out into a 10x100x100 element cube.
Now, let's look at some of the other functions available to apply to
arrays. If you remember from earlier, all of the standard arithmetic
operations are overloaded to operate on arrays one element at a time. But,
most matrix programming languages use the multiplication element to mean
matrix multiplication. This is something to keep in mind when you start
using Python. To get a true matrix multiplication, you need to use the
dot()
function.
If you have two matrices, A and B, you can multiply
them with numpy.dot(A, B)
.
Many of the standard mathematical functions,
like cosine, sine, square root and so on, are provided by NumPY and are
called universal functions. Just like with the arithmetic operators,
these universal functions are applied elementwise across the array. In
science, several common functions are used. You can get
the transpose of a matrix with the object function
a.transpose()
. If
you need to get the inverse of a matrix, there is the module function
inv()
, so you would execute
numpy.inv(a)
. The trace is also a module
function, given by numpy.trace(a)
.
Even more complicated
functions are available. You can solve systems of equations with NumPY. If
you have a matrix of coefficients given by a, and a vector of numbers
representing the righthand side of your equations given by y, you can
solve this system with numpy.solve(a,y)
. In many
situations, you may
be interested in finding the eigenvalues and eigenfunctions of a given
system. If so, you can use numpy.eig(array1)
to get
those values.
The last thing I want to look at here is a class that provides even more shortcuts, at the cost of more restrictions. Matrices (2D arrays) are so prevalent that NumPY provides a special class to optimize operations using them as much as possible. To create a new matrix, say a 2x2 matrix, you would write:
A = numpy.matrix('1.0, 2.0; 3.0, 4.0')
Now, you can get the transpose with just A.T
. Similarly, the inverse
is found with A.I
. The multiplication operation will do actual matrix
multiplication when you use matrix objects. So, given two matrix objects
A and B, you can do matrix multiplication with A*B
. The solve function
still works as expected on systems of equations that are defined using
matrix objects. Lots of tips and tricks available on the NumPY
Web site, which is well worth a look, especially as you start out.
This short introduction should get you started in thinking of Python as a viable possibility in "real" numerical computations. The NumPY module provides a very strong foundation to build up complex scientific workflows. Next month, I'll look at one of the available modules, SciPY. Until then, play with all of the raw numbercrunching possibilities provided by NumPY.
Joey Bernard has a background in both physics and computer science. This serves him well in his day job as a computational research consultant at the University of New Brunswick. He also teaches computational physics and parallel programming.