SciPy - Library of scientific algorithms for Python

Modified from J.R. Johansson ([email protected]) http://dml.riken.jp/~rob/

In [1]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib

Introduction

The SciPy framework builds on top of the low-level NumPy framework for multidimensional arrays, and provides a large number of higher-level scientific algorithms. Some of the topics that SciPy covers are:

Each of these submodules provides a number of functions and classes that can be used to solve problems in their respective topics.

In this lecture we will look at how to use some of these subpackages.

To access the SciPy package in a Python program, we start by importing everything from the scipy module.

In [2]:
from scipy import *

If we only need to use part of the SciPy framework we can selectively include only those modules we are interested in. For example, to include the linear algebra package under the name la, we can do:

In [3]:
import scipy.linalg as la

Fourier transform

Fourier transforms are one of the universal tools in computational physics, which appear over and over again in different contexts. SciPy provides functions for accessing the classic FFTPACK library from NetLib, which is an efficient and well tested FFT library written in FORTRAN. The SciPy API has a few additional convenience functions, but overall the API is closely related to the original FORTRAN library.

To use the fftpack module in a python program, include it using:

In [4]:
from scipy.fftpack import *

To demonstrate how to do a fast Fourier transform with SciPy, let's look at the following example:

In [5]:
fs = 8000.0 # sampling frequency
t = arange(1 * fs) # 1 second duration
x = sin(t / fs * 2 * pi * 1000) # 1000 Hz signal
(fig, ax) = subplots(figsize=(20, 3))
ax.plot(x[:fs * 0.010]) # look at the first 10ms
Out[5]:
[<matplotlib.lines.Line2D at 0x113c68ad0>]
In [6]:
N = 512
# calculate the fast fourier transform
X = fft(x, N)
# calculate the frequencies for the components in F
f = fftfreq(N, 1. / fs)
(fig, ax) = subplots(figsize=(20, 3))
ax.plot(f, abs(X)); # magnitude spectrum
xlabel('frequency (Hz)'); ylabel('magnitude')
Out[6]:
<matplotlib.text.Text at 0x114d516d0>

Since the signal is real, the spectrum is symmetric. We therefore only need to calculate and plot the part that corresponds to the positive frequencies:

In [7]:
X = rfft(x, N)
(fig, ax) = subplots(figsize=(10, 3))
ax.plot(linspace(0, fs / 2, N) , log(abs(X)))
xlabel('frequency (Hz)'); ylabel('log-magnitude')
Out[7]:
<matplotlib.text.Text at 0x114d750d0>

As expected, we see a peak in the spectrum centered around 1000 Hz.

Linear algebra

The linear algebra module contains a lot of matrix related functions, including linear equation solving, eigenvalue solvers, matrix functions (for example matrix-exponentiation), a number of different decompositions (SVD, LU, cholesky), etc.

Detailed documetation is available at: http://docs.scipy.org/doc/scipy/reference/linalg.html

Here we will look at how to use some of these functions:

Linear equation systems

Linear equation systems on the matrix form

$A x = b$

where $A$ is a matrix and $x,b$ are vectors can be solved like:

In [8]:
A = array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
b = array([1,2,3])
In [9]:
x = solve(A, b)
x
Out[9]:
array([-0.33333333,  0.66666667,  0.        ])
In [10]:
# check
dot(A, x)
Out[10]:
array([ 1.,  2.,  3.])

We can also do the same with

$A X = B$

where $A, B, X$ are matrices:

In [11]:
A = rand(3, 3)
B = rand(3, 3)
In [12]:
X = solve(A, B)
In [13]:
# check
norm(dot(A, X) - B)
Out[13]:
8.7419048378076915e-16

Eigenvalues and eigenvectors

The eigenvalue problem for a matrix $A$:

$\displaystyle A v_n = \lambda_n v_n$

where $v_n$ is the $n$th eigenvector and $\lambda_n$ is the $n$th eigenvalue.

To calculate eigenvalues of a matrix, use eigvals, and for calculating both eigenvalues and eigenvectors, use the function eig:

In [14]:
A = array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
evals = eigvals(A)
evals
Out[14]:
array([  1.61168440e+01,  -1.11684397e+00,  -1.30367773e-15])
In [15]:
(evals, evecs) = eig(A)
(evals, evecs)
Out[15]:
(array([  1.61168440e+01,  -1.11684397e+00,  -1.30367773e-15]),
 array([[-0.23197069, -0.78583024,  0.40824829],
        [-0.52532209, -0.08675134, -0.81649658],
        [-0.8186735 ,  0.61232756,  0.40824829]]))

The eigenvectors corresponding to the $n$th eigenvalue (stored in evals[n]) is the $n$th column in evecs, i.e., evecs[:,n]. To verify this, let's try multiplying eigenvectors with the matrix and compare to the product of the eigenvector and the eigenvalue:

In [16]:
n = 0
dot(A, evecs[:, n]), evals[n] * evecs[:, n]
Out[16]:
(array([ -3.73863537,  -8.46653421, -13.19443305]),
 array([ -3.73863537,  -8.46653421, -13.19443305]))

There are also more specialized eigensolvers, like the eigh for Hermitian matrices.

Sparse matrices

Sparse matrices are often useful in numerical simulations dealing with large systems, if the problem can be described in matrix form where the matrices or vectors mostly contains zeros. Scipy has a good support for sparse matrices, with basic linear algebra operations (such as equation solving, eigenvalue calculations, etc).

There are many possible strategies for storing sparse matrices in an efficient way. Some of the most common are the so-called coordinate form (COO), list of list (LIL) form, and compressed-sparse column CSC (and row, CSR). Each format has some advantanges and disadvantages. Most computational algorithms (equation solving, matrix-matrix multiplication, etc) can be efficiently implemented using CSR or CSC formats, but they are not so intuitive and not so easy to initialize. So often a sparse matrix is initially created in COO or LIL format (where we can efficiently add elements to the sparse matrix data), and then converted to CSC or CSR before used in real calcalations.

For more information about these sparse formats, see e.g. http://en.wikipedia.org/wiki/Sparse_matrix

When we create a sparse matrix we have to choose which format it should be stored in. For example,

In [17]:
from scipy.sparse import *
In [18]:
# dense matrix
M = array([[1, 0, 0, 0], [0, 3, 0, 0], [0, 1, 1, 0], [1, 0, 0, 1]])
M
Out[18]:
array([[1, 0, 0, 0],
       [0, 3, 0, 0],
       [0, 1, 1, 0],
       [1, 0, 0, 1]])
In [19]:
# convert from dense to sparse
A = csr_matrix(M)
A
Out[19]:
<4x4 sparse matrix of type '<type 'numpy.int64'>'
	with 6 stored elements in Compressed Sparse Row format>
In [20]:
# convert from sparse to dense
A.todense()
Out[20]:
matrix([[1, 0, 0, 0],
        [0, 3, 0, 0],
        [0, 1, 1, 0],
        [1, 0, 0, 1]])

More efficient way to create sparse matrices: create an empty matrix and populate with using matrix indexing (avoids creating a potentially large dense matrix)

In [21]:
A = lil_matrix((4, 4)) # empty 4x4 sparse matrix
A[0, 0] = 1
A[1, 1] = 3
A[2, 2] = A[2, 1] = 1
A[3, 3] = A[3, 0] = 1
A
Out[21]:
<4x4 sparse matrix of type '<type 'numpy.float64'>'
	with 6 stored elements in LInked List format>
In [22]:
A.todense()
Out[22]:
matrix([[ 1.,  0.,  0.,  0.],
        [ 0.,  3.,  0.,  0.],
        [ 0.,  1.,  1.,  0.],
        [ 1.,  0.,  0.,  1.]])

Converting between different sparse matrix formats:

In [23]:
A
Out[23]:
<4x4 sparse matrix of type '<type 'numpy.float64'>'
	with 6 stored elements in LInked List format>
In [24]:
A = csr_matrix(A)
A
Out[24]:
<4x4 sparse matrix of type '<type 'numpy.float64'>'
	with 6 stored elements in Compressed Sparse Row format>
In [25]:
A = csc_matrix(A)
A
Out[25]:
<4x4 sparse matrix of type '<type 'numpy.float64'>'
	with 6 stored elements in Compressed Sparse Column format>

We can compute with sparse matrices like with dense matrices:

In [26]:
A.todense()
Out[26]:
matrix([[ 1.,  0.,  0.,  0.],
        [ 0.,  3.,  0.,  0.],
        [ 0.,  1.,  1.,  0.],
        [ 1.,  0.,  0.,  1.]])
In [27]:
(A * A).todense()
Out[27]:
matrix([[ 1.,  0.,  0.,  0.],
        [ 0.,  9.,  0.,  0.],
        [ 0.,  4.,  1.,  0.],
        [ 2.,  0.,  0.,  1.]])
In [28]:
(A * A.T).todense()
Out[28]:
matrix([[ 1.,  0.,  0.,  1.],
        [ 0.,  9.,  3.,  0.],
        [ 0.,  3.,  2.,  0.],
        [ 1.,  0.,  0.,  2.]])
In [29]:
v = array([1, 2, 3, 4])[:, newaxis]
v
Out[29]:
array([[1],
       [2],
       [3],
       [4]])
In [30]:
# sparse matrix - dense vector multiplication
A * v
Out[30]:
array([[ 1.],
       [ 6.],
       [ 5.],
       [ 5.]])
In [31]:
# same result with dense matrix - dense vector multiplcation
A.todense() * v
Out[31]:
matrix([[ 1.],
        [ 6.],
        [ 5.],
        [ 5.]])

Optimization

Optimization (finding minima or maxima of a function) is a large field in mathematics, and optimization of complicated functions or in many variables can be rather involved. Here we will only look at a few very simple cases. For a more detailed introduction to optimization with SciPy see: http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html

To use the optimization module in scipy first include the optimize module:

In [32]:
from scipy import optimize

Finding a minima

Let's first look at how to find the minima of a simple function of a single variable:

In [33]:
def f(x):
    return 4 * x ** 3 + (x - 2) ** 2 + x ** 4
In [34]:
(fig, ax) = subplots()
x = linspace(-5, 3, 100)
ax.plot(x, f(x));

We can use the fmin_bfgs function to find the local minima of a function:

In [35]:
optimize.fmin_bfgs(f, -2)
Optimization terminated successfully.
         Current function value: -3.506641
         Iterations: 6
         Function evaluations: 30
         Gradient evaluations: 10
Out[35]:
array([-2.67298164])
In [36]:
optimize.fmin_bfgs(f, 0.5)
Optimization terminated successfully.
         Current function value: 2.804988
         Iterations: 3
         Function evaluations: 15
         Gradient evaluations: 5
Out[36]:
array([ 0.46961745])

We can also use the brent or fminbound functions. They have a bit different syntax and use different algorithms.

In [37]:
optimize.brent(f)
Out[37]:
0.46961743402759754
In [38]:
optimize.fminbound(f, -4, 2)
Out[38]:
-2.6729822917513886

Finding a solution to a function

To find the root for a function of the form $f(x) = 0$ we can use the fsolve function. It requires an initial guess:

In [39]:
omega_c = 3.0

def f(omega):
    # a transcendental equation: resonance frequencies of a low-Q 
    # SQUID terminated microwave resonator
    return tan(2 * pi * omega) - omega_c / omega
In [40]:
(fig, ax) = subplots(figsize=(10,4))
x = linspace(0, 3, 1000)
y = f(x)
mask = where(abs(y) > 50)
x[mask] = y[mask] = NaN # get rid of vertical line during plotting when the function flips sign
ax.plot(x, y)
ax.plot([0, 3], [0, 0], 'k')
ax.set_ylim(-5, 5);
-c:6: RuntimeWarning: divide by zero encountered in divide
In [41]:
optimize.fsolve(f, 0.1)
Out[41]:
array([ 0.23743014])
In [42]:
optimize.fsolve(f, 0.6)
Out[42]:
array([ 0.71286972])
In [43]:
optimize.fsolve(f, 1.1)
Out[43]:
array([ 1.18990285])

Interpolation

Interpolation is simple and convenient in scipy: The interp1d function, when given arrays describing X and Y data, returns an object that behaves like a function that can be called for an arbitrary value of x (in the range covered by X), and it returns the corresponding interpolated y value:

In [44]:
from scipy.interpolate import *
In [45]:
def f(x):
    return sin(x)
In [46]:
n = arange(0, 10)  
x = linspace(0, 9, 100)

y_meas = f(n) + 0.2 * randn(len(n)) # simulate measurement with noise
y_real = f(x)

linear_interpolation = interp1d(n, y_meas)
y_interp1 = linear_interpolation(x)

cubic_interpolation = interp1d(n, y_meas, kind='cubic')
y_interp2 = cubic_interpolation(x)
In [47]:
(fig, ax) = subplots(figsize=(20,4))
ax.plot(n, y_meas, 'bs', label='noisy data')
ax.plot(x, y_real, 'k', lw=2, label='true function')
ax.plot(x, y_interp1, 'r', label='linear interp')
ax.plot(x, y_interp2, 'g', label='cubic interp')
ax.legend(loc=3);

Further reading