Performance Python¶
Dr. Yves J. Hilpisch
The Python Quants GmbH
PyData, New York City – 08. November 2013
Performance Python¶
Python for Performance Computing:
- measure: getting to results faster
- pillar 1: agile interactive analytics & prototyping
- pillar 2: rapid application development & deployment
- pillar 3: performant execution of algorithms & code
This tutorial focuses on
- some fundamental Python issues when it comes to performance
- selected approaches to improve the performance of algorithms & code
- simple, interactive examples
It does not address such important issues like
- architectural issues regarding hardware and software
- development processes, testing, documentation and production
- real world problem modeling
Some Python Fundamentals¶
Fundamental Python Libraries¶
A fundamental Python stack for interactive data analytics and visualization should at least contain the following libraries tools:
- Python – the Python interpreter itself
- NumPy – high performance, flexible array structures and operations
- SciPy – collection of scientific modules and functions (e.g. for regression, optimization, integration)
- pandas – time series and panel data analysis and I/O
- PyTables – hierarchical, high performance database (e.g. for out-of-memory analytics)
- matplotlib – 2d and 3d visualization
- IPython – interactive data analytics, visualization, publishing
It is best to use either the Python distribution Anaconda or the Web-based analytics environment Wakari. Both provide almost "complete" Python environments.
When it comes to performance critical operations, a number of additional libraries can be used to speed-up code execution.
- numexpr – for fast numerical operations
- Numba – for dynamically compiling Python code for the CPU
- NumbaPro – for dynamically compiling Python code for multi-core CPU and GPUs
- Cython – for merging Python with C paradigms for static compilation
- IPython.Parallel – for the parallel execution of code/functions over a cluster
- ... many others
All these are also part of Anaconda and Wakari.
Comparing Performance of Python Functions¶
We define a convenience function to compare the performance of a set of Python functions executed on the same data or on different data.
def perf_comp_data(func_list, data_list, rep=3, number=1):
from timeit import repeat
res_list = {}
for name in enumerate(func_list):
stmt = name[1] + '(' + data_list[name[0]] + ')'
setup = "from __main__ import " + name[1] + ', ' + data_list[name[0]]
results = repeat(stmt=stmt, setup=setup, repeat=rep, number=number)
res_list[name[1]] = sum(results) / rep
res_sort = sorted(res_list.iteritems(), key=lambda (k, v): (v, k))
for item in res_sort:
rel = item[1] / res_sort[0][1]
print 'function: ' + item[0] + ', av. time sec: %9.5f, ' % item[1] \
+ 'relative: %6.1f' % rel
Speeding Up Numerical Computations¶
A typical compute intensive operation is to evaluate a complex mathematical expression on a large array of data.
First, let us define an example function/numerical expression.
from math import *
def f(x):
return abs(cos(x)) ** 0.5 + sin(2 + 3 * x)
Second, as our benchmark case a pure Python implementation.
I = 200000
a_py = range(I)
def f1(a):
res = []
for x in a:
res.append(f(x))
return res
One can also use different paradigms to implement the same function, like iterators or the eval function.
def f2(a):
return [f(x) for x in a]
def f3(a):
ex = 'abs(cos(x)) ** 0.5 + sin(2 + 3 * x)'
return [eval(ex) for x in a]
Third, a vectorized implementations with NumPy.
import numpy as np
a_np = np.arange(I)
def f4(a):
return (np.abs(np.cos(a)) ** 0.5 +
np.sin(2 + 3 * a))
Fourth, we use numexpr to evaluate the numerical expression.
import numexpr as ne
def f5(a):
ex = 'abs(cos(a)) ** 0.5 + sin(2 + 3 * a)'
ne.set_num_threads(1)
return ne.evaluate(ex)
def f6(a):
ex = 'abs(cos(a)) ** 0.5 + sin(2 + 3 * a)'
ne.set_num_threads(4)
return ne.evaluate(ex)
Now, let's compare the performance of all different function implementations.
func_list = ['f1', 'f2', 'f3', 'f4', 'f5', 'f6']
data_list = ['a_py', 'a_py', 'a_py', 'a_np', 'a_np', 'a_np']
perf_comp_data(func_list, data_list, rep=1)
Speeding Up Code Through Memory Layout¶
Depending on the respective algorithm you implement with NumPy, memory layout can play a non-negligible role.
import numpy as np
np.zeros((3, 3), dtype=np.float64, order='C')
The way you initialize NumPy array can have a significant influence on the performance of operations on these arrays (given a certain size of the array).
- shape: either an int, a sequence of ints or a reference to another numpy.ndarray
- dtype (optional): a numpy.dtype which are NumPy-specific basic data types for numpy.ndarray objects
- order (optional): the order of storing elements in memory, 'C' for C-like (i.e. row-wise) or 'F' for Fortran-like (i.e. column-wise)
Consider the C-like, i.e. row-wise, storage.
c = np.array([[ 1., 1., 1.],
[ 2., 2., 2.],
[ 3., 3., 3.]], order='C')
In this case, the 1s, the 2s and the 3s are stored next to each other.
By contrast, consider the Fortran-like, i.e. column-wise, storage.
f = np.array([[ 1., 1., 1.],
[ 2., 2., 2.],
[ 3., 3., 3.]], order='F')
Now, the data is stored in a way that 1, 2, 3 are next to each other for each column.
Let's see whether the memory layout makes a difference in some way when the array is large.
x = np.random.standard_normal((3, 1500000))
C = np.array(x, order='C')
F = np.array(x, order='F')
x = 0.0
Let's implement standard operations on the C-like layout array.
First, calculating sums.
%timeit C.sum(axis=0)
%timeit C.sum(axis=1)
Second, standard deviations.
%timeit C.std(axis=0)
%timeit C.std(axis=1)
Now the Fortran-like layout. Sums first.
%timeit F.sum(axis=0)
%timeit F.sum(axis=1)
Standard deviations second.
%timeit F.std(axis=0)
%timeit F.std(axis=1)
Speeding Up I/O Operations¶
SQL Databases¶
Currently, and in the foreseeable future as well, most corporate data is stored in SQL databases. We take this as a starting point and generate a SQL table with dummy data.
filename = 'numbers'
import sqlite3 as sq3 # imports the SQLite library
# Remove files if they exist
import os
try:
os.remove(filename + '.db')
os.remove(filename + '.h5')
os.remove(filename + '.hs')
os.remove(filename + '.csv')
except:
pass
We create a table and populate the table with the dummy data.
query = 'CREATE TABLE numbs (no1 real, no2 real, no3 real, no4 real, no5 real)'
con = sq3.connect(filename + '.db')
con.execute(query)
con.commit()
# Writing Random Numbers to SQLite3 Database
rows = 5000000 # 5,000,000 rows for the SQL table (ca. 250 MB in size)
nr = np.random.standard_normal((rows, 5))
%time con.executemany('INSERT INTO numbs VALUES(?, ?, ?, ?, ?)', nr)
con.commit()
We can now read the data into a NumPy array.
# Reading a couple of rows from SQLite3 database table
array(con.execute('SELECT * FROM numbs').fetchmany(3)).round(3)
# Reading the whole table into an NumPy array
%time result = array(con.execute('SELECT * FROM numbs').fetchall())
PyTables, which is based on the HDF5 file format, allows much faster, disk-based I/O than a typical SQL database. For example, writing the same data is much faster compared to SQLite3.
import tables as tb
%%time
h5 = tb.open_file('numbers.h5', 'w')
h5.create_array('/', 'numbs', obj=result, title='numbs')
h5.close()
Reading speed is equally impressive.
%%time
h5 = tb.open_file('numbers.h5', 'r')
arr = h5.root.numbs.read()
h5.close()
arr.nbytes
Reading and Writing Data with pandas¶
pandas is optimized both for convenience and performance when it comes to analytics and I/O operations.
import pandas as pd
import pandas.io.sql as pds
%time data = pds.read_frame('SELECT * FROM numbs', con)
con.close()
data.head()
The plotting of larger data sets often is a problem. This can easily be accomplished here by only plotting every 500th data row, for instance.
%time data[::500].cumsum().plot(grid=True, legend=True) # every 500th row is plotted
There are some reasons why you may want to write the DataFrame to disk:
- Database Performance: you do not want to query the SQL database everytime anew
- IO Speed: you want to increase data/information retrievel
- Data Structure: maybe you want to store additional data with the original data set or you have generated a DataFrame out of multiple sources
pandas is tightly integrated with PyTables. It makes data storage and I/O operations with pandas not only convenient but quite fast as well.
%%time
h5 = pd.HDFStore(filename + '.hs', 'w')
h5['data'] = data # writes the DataFrame to HDF5 database
h5.close()
Now, we can read the data from the HDF5 file into another DataFrame object. Performance of such an operation generally is quite high.
%%time
h5 = pd.HDFStore(filename + '.hs', 'r')
data_new = h5['data']
h5.close()
%time data.sum() - data_new.sum() # checking if the DataFrames are the same
Some final checks.
ll numb*
Speeding Up Code Execution through Distributed Computing¶
The algorithm we (arbitrarily) choose is the Monte Carlo estimator for a European option value in the Black-Scholes-Merton set-up. In this set-up the underlying follows the stochastic differential equation:
\[dS_t = r S_t dt + \sigma S_t dZ_t\]
The Monte Carlo estimator for a call option value looks like follows:
\(C_0 = e^{-rT} \frac{1}{I} \sum_{i=1}^{I} \max[S_T(i) - K, 0]\)
where \(S_T(i)\) is the \(i\)-th simulated value of the underlying at maturity \(T\).
A dynamic implementation of the BSM Monte Carlo valuation could look like follows.
def BSM_MCS(S0):
''' Dynamic Black-Scholes-Merton MCS Estimator for European Calls. '''
import numpy as np
K = 100.; T = 1.0; r = 0.05; vola = 0.2
M = 50; I = 20000
dt = T / M
rand = np.random.standard_normal((M + 1, I))
S = np.zeros((M + 1, I)); S[0] = S0
for t in range(1, M + 1):
S[t] = S[t-1] * np.exp((r - 0.5 * vola ** 2) * dt + vola * np.sqrt(dt) * rand[t])
BSM_MCS_Value = np.sum(np.maximum(S[-1] - K, 0)) / I
return BSM_MCS_Value
The Sequential Calculation¶
The benchmark case of sequential valuation of \(n\) options.
def seqValue(n):
optionValues = []
for S in linspace(80, 100, n):
optionValues.append(BSM_MCS(S))
return optionValues
The Parallel Calculation¶
First, start a (local) cluster.
# in the shell ...
####
# ipcluster start -n 4
####
# or maybe more cores
Second, enable parallel computing capabilities.
from IPython.parallel import Client
c = Client(profile="default")
view = c.load_balanced_view()
The function for the parallel execution of a number \(n\) of different options valuations.
def parValue(n):
optionValues = []
for S in linspace(80, 100, n):
value = view.apply_async(BSM_MCS, S)
# asynchronously calculate the option values
optionValues.append(value)
c.wait(optionValues)
return optionValues
Let us compare performance of the two different approaches.
n = 20 # number of option valuations
func_list = ['seqValue', 'parValue']
data_list = 2 * ['n']
perf_comp_data(func_list, data_list, rep=3)
Performance scales linearly with the number of CPU cores available.
Speeding Up Code through Dynamic Compiling¶
Numba is an open-source, NumPy-aware optimizing compiler for Python sponsored by The Python Quants GmbH It uses the remarkable LLVM compiler infrastructure to compile Python byte-code to machine code especially for use in the NumPy run-time and SciPy modules.
Introductory Examples¶
Let us start with a problem that typically lead to performance issues in Python: alogrithms with a nested loops. A sandbox variant can illustrate the problem.
from math import cos, log
def f_py(I, J):
res = 0
for i in range(I):
for j in range (J):
res += int(cos(log(1)))
return res
I, J = 1000, 1000
%time f_py(I, J)
In principle, this can be vectorized with NumPy arrays.
def f_np(I, J):
a = np.ones((I, J), dtype=np.float64)
return int(np.sum(np.cos(np.log(a)))), a
%time res, a = f_np(I, J)
This is much faster, but not really memory efficient ...
a.nbytes
Numba can provide help.
import numba as nb
%time f_nb = nb.autojit(f_py)
It preserves memory efficiency while speeding up execution considerably in this case.
%time f_nb(I, J)
Let us compare the performance of the different alternatives.
func_list = ['f_py', 'f_np', 'f_nb']
data_list = 3 * ['I, J']
perf_comp_data(func_list, data_list, rep=3)
N-Body Simulation¶
The simulation of the movement of N bodies in relation to each other, where \(N\) is a large number and a body can be a star/planet or a small particle, is an important problem in Physics.
The basic simulation algorithm involves a number of iterations ("loops") since one has to determine the force that every single body exerts on all the other bodies in the system to be simulated (e.g. a galaxy or a system of atomic particles).
The algorithm lends iteself pretty well for parallization. However, the example code that follows does not use any parallelization techniques. It illustrates rather what speed-ups are achievable by
- using vectorization approaches with NumPy and
- just-in-time compiling capabilities of Numba
A description and discussion of the algorithm is found in the recent article Vladimirov, Andrey and Vadim Karpusenko (2013): "Test-Driving Intel Xeon-Phi Coprocessors with a Basic N-Body Simulation." White Paper, Colfax International.
Pure Python¶
We first use pure Python to implement the algorithm. We start by generating random vectors of particle locations and a zero vector for the speed of each particle.
import random
nParticles = 150
nSteps = 15
p_py = []
for i in range(nParticles):
p_py.append([random.gauss(0.0, 1.0) for j in range(3)])
pv_py = [[0, 0, 0] for x in p_py]
Next, we define the simulation function for the movements of each particle given their gravitational interaction.
def nbody_py(particle, particlev): # lists as input
dt = 0.01
for step in range(1, nSteps + 1, 1):
for i in range(nParticles):
Fx = 0.0; Fy = 0.0; Fz = 0.0
for j in range(nParticles):
if j != i:
dx = particle[j][0] - particle[i][0]
dy = particle[j][1] - particle[i][1]
dz = particle[j][2] - particle[i][2]
drSquared = dx * dx + dy * dy + dz * dz
drPowerN32 = 1.0 / (drSquared + sqrt(drSquared))
Fx += dx * drPowerN32
Fy += dy * drPowerN32
Fz += dz * drPowerN32
particlev[i][0] += dt * Fx
particlev[i][1] += dt * Fy
particlev[i][2] += dt * Fz
for i in range(nParticles):
particle[i][0] += particlev[i][0] * dt
particle[i][1] += particlev[i][1] * dt
particle[i][2] += particlev[i][2] * dt
return particle, particlev
We then execute the function to measure execution speed.
%time p_py, pv_py = nbody_py(p_py, pv_py)
NumPy Solution¶
The basic algorithm allows for a number of vectorizations by using NumPy.
import numpy as np
p_np = np.random.standard_normal((nParticles, 3))
pv_np = np.zeros_like(p_np)
Using NumPy arrays and matplotlib's 3d plotting capabilities, we can easily visualize the intial positions of the particles in the space.
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(p_np[:, 0], p_np[:, 1], p_np[:, 2],
c='r', marker='.')
The respective NumPy simulation function is considerably more compact.
def nbody_np(particle, particlev): # NumPy arrays as input
nSteps = 15; dt = 0.01
for step in range(1, nSteps + 1, 1):
Fp = np.zeros((nParticles, 3))
for i in range(nParticles):
dp = particle - particle[i]
drSquared = np.sum(dp ** 2, axis=1)
h = drSquared + np.sqrt(drSquared)
drPowerN32 = 1. / np.maximum(h, 1E-10)
Fp += -(dp.T * drPowerN32).T
particlev += dt * Fp
particle += particlev * dt
return particle, particlev
It is also considerably faster.
%time p_np, pv_np = nbody_np(p_np, pv_np)
Let's check how the particles moved. The following is a snapshot of the final positions of all particles.
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(p_np[:, 0], p_np[:, 1], p_np[:, 2],
c='r', marker='.')
Dynamically Compiled Version¶
For this case, we first implement a function which is somehow a hybrid function of the pure Python (loop-heavy) and the NumPy array class-based version.
def nbody_hy(particle, particlev): # NumPy arrays as input
dt = 0.01
for step in range(1, nSteps + 1, 1):
for i in range(nParticles):
Fx = 0.0; Fy = 0.0; Fz = 0.0
for j in range(nParticles):
if j != i:
dx = particle[j,0] - particle[i,0]
dy = particle[j,1] - particle[i,1]
dz = particle[j,2] - particle[i,2]
drSquared = dx * dx + dy * dy + dz * dz
drPowerN32 = 1.0 / (drSquared + sqrt(drSquared))
Fx += dx * drPowerN32
Fy += dy * drPowerN32
Fz += dz * drPowerN32
particlev[i, 0] += dt * Fx
particlev[i, 1] += dt * Fy
particlev[i, 2] += dt * Fz
for i in range(nParticles):
particle[i,0] += particlev[i,0] * dt
particle[i,1] += particlev[i,1] * dt
particle[i,2] += particlev[i,2] * dt
return particle, particlev
This function is then compiled with Numba.
import numba as nb
nbody_nb = nb.autojit(nbody_hy)
This new, compiled version of the hybrid function is yet faster.
p_nb = np.random.standard_normal((nParticles, 3))
pv_nb = np.zeros_like(p_nb)
%time p_nb, pv_nb = nbody_nb(p_nb, pv_nb)
Speed Comparisons¶
The dynamically compiled version is by far the fastest.
func_list = ['nbody_py', 'nbody_np', 'nbody_nb', 'nbody_hy']
data_list = ['p_py, pv_py', 'p_np, pv_np', 'p_nb, pv_nb', 'p_nb, pv_nb']
perf_comp_data(func_list, data_list, rep=3)
Performance Python Approaches¶
Nowadays, the Python ecosystem provides a number of ways to improve the performance of code.
- paradigms: some Python paradigms might be more performant than others given a specific problem
- libraries: there is a wealth of libraries available for different types of problems which often lead to much higher performance given a problem that fits into the scope of the library
- compiling: a number of powerful compiling solutions are available, static and dynamic ones
- parallelization: some Python libraries have built-in parallellization capabilities, others allow to harness the full power of multiple cores CPUs, whole clusters or GPUs
A major benefit of the Python ecosystem is, that all these approaches generally are easily implementable, meaning that the additional effort included is generally quite low (even for non-experts). In other words, performance improvements often are low hanging fruits.
The Python Quants GmbH ¶
The Python Quants – the company Web site
Dr. Yves J. Hilpisch – my personal Web site
Derivatives Analytics with Python – my new book
Contact Us