Solving multiple linear sparse matrix equations: numpy.linalg.solve vs. scipy.sparse.linalg.spsolve

2024/4/15 0:43:33

I have to solve a large amount of linear matrix equations of the type "Ax=B" for x where A is a sparse matrix with mainly the main diagonal populated and B is a vector.

My first approach was to use dense numpy arrays for this purpose with numpy.linalg.solve, and it works fine with a (N,n,n)-dimensional array with N being the number of linear matrix equations and n the square matrix dimension. I first used it with a for loop iterating through all equations, which in fact is rather slow. But then realized that you can also pass the (N,n,n)-dimensional matrix directly to numpy.linalg.solve without any for loop (which by the way I did not find in the documentation I read). This already gave a good increase in computation speed (details see below).

However, because I have sparse matrices, I also had a look at the scipy.sparse.linalg.spsolve function which does similar things like the corresponding numpy function. Using a for loop iterating through all equations is much, much faster than the numpy solution, but it seems impossible to pass the (N,n,n)-dimesional array directly to scipy´s spsolve.

Here is what I tried so far:

First, I calculate some fictional A matrices and B vectors with random numbers for test purposes, both sparse and dense:

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolvenumber_of_systems = 100 #corresponds to N in the text
number_of_data_points = 1000 #corresponds to n in the text#calculation of sample matrices (dense and sparse)
A_sparse = np.empty(number_of_systems,dtype=object)
A_dense = np.empty((number_of_systems,number_of_data_points,number_of_data_points))for ii in np.arange(number_of_systems):A_sparse[ii] = sparse.spdiags(np.random.random(number_of_data_points),0,number_of_data_points,number_of_data_points)A_dense[ii] = A_sparse[ii].todense()#calculation of sample vectors
B = np.random.random((number_of_systems,number_of_data_points))

1) First approach: numpy.linalg.solve with for loop:

def solve_dense_3D(A,B):results = np.empty((A.shape[0],A.shape[1]))for ii in np.arange(A.shape[0]):results[ii] = np.linalg.solve(A[ii],B[ii])return resultsresult_dense_for = solve_dense_3D(A_dense,B)

Timing:

timeit(solve_dense_3D(A_dense,B))
1.25 s ± 27.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

2) Second approach: Passing the (N,n,n)-dimensional matrix directly to numpy.linalg.solve:

result_dense = np.linalg.solve(A_dense,B)

Timing:

timeit(np.linalg.solve(A_dense,B))
769 ms ± 9.68 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

3) Third approach: Using scipy.sparse.linalg.spsolve with a for loop:

def solve_sparse_3D(A,B):results = np.empty((A.shape[0],A[0].shape[0]))for ii in np.arange(A.shape[0]):results[ii] = spsolve(A[ii],B[ii])return resultsresult_sparse_for = solve_sparse_3D(A_sparse,B)

Timing:

timeit(solve_sparse_3D(A_sparse,B))
30.9 ms ± 132 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

It is obvoius that there is an advantage being able to omit the for loop from approach 1 and 2. By far the fastest alternative is, as could probably be expected, approach 3 with sparse matrices.

Because the number of equations in this example is still rather low for me and because I have to do things like that very often, I would be happy if there was a solution using scipy´s sparse matrices without a for loop. Is anybody aware of a way to achieve that? Or maybe there is another way to solve the problem in an even different way? I would be happy for suggestions.

Answer

A small demo outlining the idea from my comment above:

""" YOUR CODE (only imports changed + deterministic randomness) """import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve
from scipy.sparse import block_diag
from time import perf_counter as pcnp.random.seed(0)number_of_systems = 100 #corresponds to N in the text
number_of_data_points = 1000 #corresponds to n in the text#calculation of sample matrices (dense and sparse)
A_sparse = np.empty(number_of_systems,dtype=object)
A_dense = np.empty((number_of_systems,number_of_data_points,number_of_data_points))for ii in np.arange(number_of_systems):A_sparse[ii] = sparse.spdiags(np.random.random(number_of_data_points),0,number_of_data_points,number_of_data_points)A_dense[ii] = A_sparse[ii].todense()#calculation of sample vectors
B = np.random.random((number_of_systems,number_of_data_points))def solve_sparse_3D(A,B):results = np.empty((A.shape[0],A[0].shape[0]))for ii in np.arange(A.shape[0]):results[ii] = spsolve(A[ii],B[ii])return resultsstart = pc()
result_sparse_for = solve_sparse_3D(A_sparse,B)
end = pc()
print(result_sparse_for)
print(end - start)""" ALTERNATIVE APPROACH """def solve_sparse_3D_blockdiag(A,B):oldN = B.shape[0]A_ = block_diag(A)    # huge sparse block-matrix of independent problemsB_ = B.ravel()        # flattened vectorresults = spsolve(A_, B_)return results.reshape(oldN, -1)    # unflatten resultsstart = pc()
result_sparse_for = solve_sparse_3D_blockdiag(A_sparse,B)
end = pc()
print(result_sparse_for)
print(end - start)

which outputs

[[ 0.97529866  1.26406276  0.83348888 ...  0.99310639  3.90781207
0.16724226]
[ 1.23398934 28.82088739  1.6955886  ...  1.85011686  0.23386882
1.17208753]
[ 0.92864777  0.22248781  0.09445412 ...  2.5080376   0.91701228
0.97266564]
...
[ 0.33087093  0.89034736  1.7523883  ...  0.2171746   4.89236164
0.31546549]
[ 1.2163625   3.0100941   0.87216264 ...  1.62105596  0.33211353
2.07929302]
[ 5.35677404  1.23830776  0.16073721 ...  0.26492506  0.53676822
3.73192617]]
0.08764066299999995###[[ 0.97529866  1.26406276  0.83348888 ...  0.99310639  3.90781207
0.16724226]
[ 1.23398934 28.82088739  1.6955886  ...  1.85011686  0.23386882
1.17208753]
[ 0.92864777  0.22248781  0.09445412 ...  2.5080376   0.91701228
0.97266564]
...
[ 0.33087093  0.89034736  1.7523883  ...  0.2171746   4.89236164
0.31546549]
[ 1.2163625   3.0100941   0.87216264 ...  1.62105596  0.33211353
2.07929302]
[ 5.35677404  1.23830776  0.16073721 ...  0.26492506  0.53676822
3.73192617]]
0.07241856000000013

There are some things to do:

  • use your original more sane benchmarking-approach
  • build the block_diag in the correct sparse format to get rid of some potential warning and slowdown: see docs
  • tune spsolve's parameter permc_spec
https://en.xdnf.cn/q/73270.html

Related Q&A

I want to return html in a flask route [duplicate]

This question already has answers here:Python Flask Render Text from Variable like render_template(4 answers)Closed 6 years ago.Instead of using send_static_file, I want to use something like html(<…

Why doesnt cv2 dilate actually affect my image?

So, Im generating a binary (well, really gray scale, 8bit, used as binary) image with python and opencv2, writing a small number of polygons to the image, and then dilating the image using a kernel. Ho…

How to plot text clusters?

I have started to learn clustering with Python and sklearn library. I have wrote a simple code for clustering text data. My goal is to find groups / clusters of similar sentences. I have tried to plot…

Selenium - Unresponsive Script Error (Firefox)

This question has been asked before, but the answer given does not seem to work for me. The problem is, when opening a page using Selenium, I get numerous "Unresponsive Script" pop ups, refe…

Fail to validate URL in Facebook webhook subscription with python flask on the back end and ssl

Im trying to start using new messenger platform from FB. So i have server with name (i.e.) www.mysite.com I got a valid SSL certificate for that domain and apache is setup correctly - all good.I have …

What is a proper way to test SQLAlchemy code that throw IntegrityError?

I have read this Q&A, and already try to catch exception on my code that raise an IntegrityError exception, this way :self.assertRaises(IntegrityError, db.session.commit())But somehow my unit test …

What are the different options for social authentication on Appengine - how do they compare?

[This question is intended as a means to both capture my findings and sanity check them - Ill put up my answer toute suite and see what other answers and comments appear.]I spent a little time trying t…

Is there any way to get source code inside context manager as string?

Source code of function can be received with inspect.getsourcelines(func) function. Is there any way to do same for context manager?with test():print(123)# How to get "print(123)" as line he…

Create temporary file in Python that will be deleted automatically after sometime

Is it possible to create temporary files in python that will be deleted after some time? I checked tempfile library which generates temporary files and directories.tempfile.TemporaryFile : This functi…

Python GTK+ Canvas

Im currently learning GTK+ via PyGobject and need something like a canvas. I already searched the docs and found two widgets that seem likely to do the job: GtkDrawingArea and GtkLayout. I need a few b…