You can solve matrices of linear systems of equations in numpy in parallel using multithreaded implementations of the algorithms.

In this tutorial, you will discover how to solve systems of linear equations using multithreaded numpy functions.

Let’s get started.

Table of Contents

## Numpy Matrix Solvers are Multithreaded

Linear algebra is a field of mathematics concerned with linear equations with arrays and matrices of numbers.

Numpy is a Python library for working with arrays of numbers. As such, it implements many linear algebra functions in the **numpy.linalg** module.

Included are functions for solving systems of linear equations.

Many of the **numpy.linalg** functions for working with matrices call down into the BLAS and LAPACK libraries. These libraries provide efficient implementations of common matrix functions. Some of these efficient implementations are multithreaded.

The NumPy linear algebra functions rely on BLAS and LAPACK to provide efficient low level implementations of standard linear algebra algorithms.

— Linear algebra (numpy.linalg)

Numpy calls down into LAPACK libraries to solve systems of linear equations. The implementation of these functions are multithreaded.

The specific solvers that are multithreaded include:

- Solve linear matrix equation:
**numpy.linalg.solve()**, which calls**_gesv**in LAPACK. - Solve linear least squares:
**numpy.linalg.lstsq()**

The BLAS and LAPACK libraries will use all CPU cores by default. Nevertheless, we can configure the libraries to use a specific number of threads, such as one thread for a single-threaded version of a function and one thread per physical CPU core, which often results in better performance.

You can learn more about BLAS and LAPACK in numpy in the tutorial:

This can be achieved using an environment variable, such as **OMP_NUM_THREADS** for OpenBLAS. Check the documentation for your BLAS library in order to discover which environment variable to set.

You can learn more about setting the number of BLAS threads in the tutorial:

Now that we know that common linear algebra matrix solvers are multithreaded, let’s look at some worked examples.

## Parallel Matrix Solver

In this section, we will explore a parallel implementation of the linear matrix equation solver.

In mathematics, a system of linear equations (or linear system) is a collection of one or more linear equations involving the same variables.

— System of linear equations, Wikipedia.

We can implement the matrix multiplication with NumPy via the **numpy.linalg.solve()** function.

Solve a linear matrix equation, or system of linear scalar equations. Computes the “exact” solution, x, of the well-determined, i.e., full rank, linear matrix equation ax = b.

— numpy.linalg.solve API

Let’s take a closer look at a single-threaded and multithreaded version of the algorithm.

### Single Threaded Matrix Solver

In this example, we will configure BLAS to use a single thread.

We will create a square matrix with the dimensions 8,000 x 8,000 and a vector with 8,000 elements to solve the system of linear equations.

The complete example is listed below.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
# example of single threaded matrix solver from os import environ environ['OMP_NUM_THREADS'] = '1' from time import time from numpy.random import rand from numpy.linalg import solve # record the start time start = time() # size of arrays n = 8000 # create matrix a = rand(n, n) # create result b = rand(n, 1) # solve least squares equation x = solve(a, b) # calculate and report duration duration = time() - start print(f'Took {duration:.3f} seconds') |

Running the example took about 7.070 seconds on my system.

It may take more or fewer seconds, depending on the hardware in your system and the versions of Python, Numpy, Scipy, and the BLAS library installed.

1 |
Took 7.070 seconds |

Next, we can run the same example with a multithreaded algorithm.

### Multithreaded Matrix Solver

We can update the above example to use one thread per physical CPU core in the system.

I have 4 cores, therefore, OpenBLAS is configured to use 4 threads. Set the number of threads to match the number of cores in your system accordingly.

The complete example is listed below.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
# example of multithreaded matrix solver from os import environ environ['OMP_NUM_THREADS'] = '4' from time import time from numpy.random import rand from numpy.linalg import solve # record the start time start = time() # size of arrays n = 8000 # create matrix a = rand(n, n) # create result b = rand(n, 1) # solve least squares equation x = solve(a, b) # calculate and report duration duration = time() - start print(f'Took {duration:.3f} seconds') |

Running the example took about 3.098 seconds.

That is about 3.972 seconds faster or a speed-up of about 2.28x.

1 |
Took 3.098 seconds |

## Parallel Matrix Least Squares Solver

In this section, we will explore a parallel implementation of the linear least squares solver.

That is, using the least squares method to solve a system of linear equations.

Linear least squares (LLS) is the least squares approximation of linear functions to data. It is a set of formulations for solving statistical problems involved in linear regression, including variants for ordinary (unweighted), weighted, and generalized (correlated) residuals.

— Linear least squares, Wikipedia.

We can implement the matrix multiplication with NumPy via the **numpy.linalg.lstsq()** function.

Return the least-squares solution to a linear matrix equation. Computes the vector x that approximately solves the equation a @ x = b.

— numpy.linalg.lstsq API.

Let’s take a closer look at a single-threaded and multithreaded version of the algorithm.

### Single Threaded Matrix Least Squares Solver

In this example, we will configure BLAS to use a single thread.

We will create a square matrix with the dimensions 5,000 x 5,000 and a vector with 5,000 elements to solve the system of linear equations.

The complete example is listed below.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
# example of single threaded matrix least squares solution from os import environ environ['OMP_NUM_THREADS'] = '1' from time import time from numpy.random import rand from numpy.linalg import lstsq # record the start time start = time() # size of arrays n = 5000 # create matrix a = rand(n, n) # create result b = rand(n, 1) # solve least squares equation x, _, _, _ = lstsq(a, b, rcond=None) # calculate and report duration duration = time() - start print(f'Took {duration:.3f} seconds') |

Running the example took about 35.347 seconds on my system.

1 |
Took 35.347 seconds |

Next, we can run the same example with a multithreaded algorithm.

### Multithreaded Matrix Least Squares Solver

We can update the above example to use one thread per physical CPU core in the system.

I have 4 cores, therefore, OpenBLAS is configured to use 4 threads. Set the number of threads to match the number of cores in your system accordingly.

The complete example is listed below.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
# example of multithreaded matrix least squares solution from os import environ environ['OMP_NUM_THREADS'] = '4' from time import time from numpy.random import rand from numpy.linalg import lstsq # record the start time start = time() # size of arrays n = 5000 # create matrix a = rand(n, n) # create result b = rand(n, 1) # solve least squares equation x, _, _, _ = lstsq(a, b, rcond=None) # calculate and report duration duration = time() - start print(f'Took {duration:.3f} seconds') |

Running the example took about 29.909 seconds.

That is about 5.438 seconds faster or a speed-up of about 1.18x.

1 |
Took 29.909 seconds |

## Takeaways

You now know how to solve systems of linear equations using multithreaded numpy functions.

**Do you have any questions?**

Ask your questions in the comments below and I will do my best to answer.

Photo by Pascal Bernardon on Unsplash

## Leave a Reply