9. Python for Scientific Computing#

9.2.1. Why do we need them?#

One reason we use scientific libraries is because they implement routines we want to use.

• numerical integration, interpolation, linear algebra, root finding, etc.

For example, it’s usually better to use an existing routine for root finding than to write a new one from scratch.

(For standard algorithms, efficiency is maximized if the community can coordinate on a common set of implementations, written by experts and tuned by users to be as fast and robust as possible!)

But this is not the only reason that we use Python’s scientific libraries.

Another is that pure Python is not fast.

So we need libraries that are designed to accelerate execution of Python code.

They do this using two strategies:

1. using compilers that convert Python-like statements into fast machine code for individual threads of logic and

2. parallelizing tasks across multiple “workers” (e.g., CPUs, individual threads inside GPUs).

We will discuss these ideas extensively in this and the remaining lectures from this series.

9.2.2. Python’s Scientific Ecosystem#

At QuantEcon, the scientific libraries we use most often are

• NumPy

• SciPy

• Matplotlib

• JAX

• Pandas

• Numba

Here’s how they fit together:

• NumPy forms foundations by providing a basic array data type (think of vectors and matrices) and functions for acting on these arrays (e.g., matrix multiplication).

• SciPy builds on NumPy by adding numerical methods routinely used in science (interpolation, optimization, root finding, etc.).

• Matplotlib is used to generate figures, with a focus on plotting data stored in NumPy arrays.

• JAX includes array processing operations similar to NumPy, automatic differentiation, a parallelization-centric just-in-time compiler, and automated integration with hardware accelerators such as GPUs.

• Pandas provides types and functions for manipulating data.

• Numba provides a just-in-time compiler that plays well with NumPy and helps accelerate Python code.

We will discuss all of these libraries extensively in this lecture series.

9.3.1. High vs low level code#

Higher-level languages like Python are optimized for humans.

This means that the programmer can leave many details to the runtime environment

• specifying variable types

• memory allocation/deallocation

• etc.

In addition, pure Python is run by an interpreter, which executes code statement-by-statement.

This makes Python flexible, interactive, easy to write, easy to read, and relatively easy to debug.

On the other hand, the standard implementation of Python (called CPython) cannot match the speed of compiled languages such as C or Fortran.

9.3.2. Where are the bottlenecks?#

Why is this the case?

a, b = 10, 10
a + b

20

a, b = 'foo', 'bar'
a + b

'foobar'

a, b = ['foo'], ['bar']
a + b

['foo', 'bar']

#include <stdio.h>

int main(void) {
    int i;
    int sum = 0;
    for (i = 1; i <= 10; i++) {
        sum = sum + i;
    }
    printf("sum = %d\n", sum);
    return 0;
}

9.3.3. Data Access#

Another drag on speed for high-level languages is data access.

To illustrate, let’s consider the problem of summing some data — say, a collection of integers.

9.3.4. Summary#

Does the discussion above mean that we should just switch to C or Fortran for everything?

The answer is: Definitely not!

For any given program, relatively few lines are ever going to be time-critical.

Hence it is far more efficient to write most of our code in a high productivity language like Python.

Moreover, even for those lines of code that are time-critical, we can now equal or outpace binaries compiled from C or Fortran by using Python’s scientific libraries.

On that note, we emphasize that, in the last few years, accelerating code has become essentially synonymous with parallelization.

This task is best left to specialized compilers!

Certain Python libraries have outstanding capabilities for parallelizing scientific code – we’ll discuss this more as we go along.

9.4.1. Vectorization#

One method for avoiding memory traffic and type checking is array programming.

Many economists usually refer to array programming as “vectorization.”

The key idea is to send array processing operations in batch to pre-compiled and efficient native machine code.

The machine code itself is typically compiled from carefully optimized C or Fortran.

For example, when working in a high level language, the operation of inverting a large matrix can be subcontracted to efficient machine code that is pre-compiled for this purpose and supplied to users as part of a package.

The core benefits are

1. type-checking is paid per array, rather than per element, and

2. arrays containing elements with the same data type are efficient in terms of memory access.

The idea of vectorization dates back to MATLAB, which uses vectorization extensively.

9.4.2. Vectorization vs for pure Python loops#

Let’s try a quick speed comparison to illustrate how vectorization can accelerate code.

Here’s some non-vectorized code, which uses a native Python loop to generate, square and then sum a large number of random variables:

n = 1_000_000

with qe.Timer():
    y = 0      # Will accumulate and store sum
    for i in range(n):
        x = random.uniform(0, 1)
        y += x**2

0.25 seconds elapsed

The following vectorized code uses NumPy, which we’ll soon investigate in depth, to achieve the same thing.

with qe.Timer():
    x = np.random.uniform(0, 1, n)
    y = np.sum(x**2)

0.01 seconds elapsed

As you can see, the second code block runs much faster.

It breaks the loop down into three basic operations

1. draw n uniforms

2. square them

3. sum them

These are sent as batch operators to optimized machine code.

9.4.3. JIT compilers#

At best, vectorization yields fast, simple code.

However, it’s not without disadvantages.

One issue is that it can be highly memory-intensive.

This is because vectorization tends to create many intermediate arrays before producing the final calculation.

Another issue is that not all algorithms can be vectorized.

Because of these issues, most high performance computing is moving away from traditional vectorization and towards the use of just-in-time compilers.

In later lectures in this series, we will learn about how modern Python libraries exploit just-in-time compilers to generate fast, efficient, parallelized machine code.

9.5.1. Parallelization on CPUs#

Let’s review the two main kinds of CPU-based parallelization commonly used in scientific computing and discuss their pros and cons.

9.5.2. Hardware Accelerators#

While CPUs with multiple cores have become standard for parallel computing, a more dramatic shift has occurred with the rise of specialized hardware accelerators.

These accelerators are designed specifically for the kinds of highly parallel computations that arise in scientific computing, machine learning, and data science.

9.5.2.1. GPUs and TPUs#

The two most important types of hardware accelerators are

• GPUs (Graphics Processing Units) and

• TPUs (Tensor Processing Units).

GPUs were originally designed for rendering graphics, which requires performing the same operation on many pixels simultaneously.

Scientists and engineers realized that this same architecture — many simple processors working in parallel — is ideal for scientific computing tasks

TPUs are a more recent development, designed by Google specifically for machine learning workloads.

Like GPUs, TPUs excel at performing massive numbers of matrix operations in parallel.

9.5.3. Single GPUs vs GPU Servers#

There are two common ways to access GPU resources:

9.5.3.2. Multi-GPU Servers#

For larger-scale problems, servers containing multiple GPUs (often 4-8 GPUs per server) are increasingly common.

With appropriate software, computations can be distributed across multiple GPUs, either within a single server or across multiple servers.

This enables researchers to tackle problems that would be infeasible on a single GPU or CPU.

9.5.4. Summary#

GPU computing is becoming far more accessible, particularly from within Python.

Some Python scientific libraries, like JAX, now support GPU acceleration with minimal changes to existing code.

We will explore GPU computing in more detail in later lectures, applying it to a range of economic applications.