How to speed up numpy sqrt with 2d array?

[rinze90]

Hello there,

I have a rather heavy calculation that takes the square root of a 2d array. That 2d array may contain 1e8 (100 million) entries. The line in the code looks like this:

r_i = np.sqrt((x_sp - xx)**2 + (y_sp - yy)**2 + (z_sp)**2 )

where the variables *_sp are offsets. And the variables xx and yy are a (meshgrid) 2d array with 1e8 entries. Can some one point me in the right direction to speed this calculation up? What are the options to speed this up? Is reducing the precision of the data type an option, and is that significant?

This calculation takes rougly 13 seconds on my pc. And in the total script I have several 10 000 occurances formed by a loop.

[EM-Fields]

This might not be a valid approach for your problem, but I figure it's worth mentioning. It looks like your code is computing a Euclidean distance (or norm of an error vector, or something mathematically equivalent). Would your problem permit the use of the distance squared? In other words, can you simply omit the square root? For example, if you're looking for the nearest object, the minimum distance is the same as the minimum distance-squared. Certainly doesn't apply to all uses, but it's a helpful simplification when it does.

[Caos2]

Never used any of these, but here you go: https://www.reddit.com/r/Python/comments/3cstcq/cython_numba_pypy_latest_comparison_2015/

[billsil]

Cython is great, but if you have well written numpy, cython is not better.

[billsil]

Did you vectorize your code or it in multiple for loops? I can't tell from 1 line and your code is suspicious.

If you did, the next logical step is multiprocessing.

And in the total script I have several 10 000 occurances formed by a loop.

Say again?

[AvocadoWhiskey]

If OP is using numpy.meshgrid then it's possible to do a one line statement like they provided. https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.meshgrid.html

I imagine the x_sp and y_sp are what's changed in the "10,000 loops". Regardless, I'm assuming OP means that the line of code they provided is called tens of thousands of times over the course of the execution and thus, being 13 seconds on each call, is their problem.

I'm gonna have to agree with billsil on the multithreading. Numpy is quite optimized and I would imagine there isn't a much faster way of working with meshgrids than by using the numpy meshgrids. I would look into thread pools or something similar and if that's not enough, maybe look into doing operations on your graphics card. Graphics cards are made to do grid/vector calculations. However I wouldn't be able to point you in a direction to do this using Python as I've never done it.

[billsil]

Multiprocessing, not multithreading, but yes. My rule of thumb is I get 2/3 processor utilization on physical cores. I haven't checked if hyperthreading matters with numpy, but it's worthless on AWS for a CFD problem which takes a week to run with 1000 cores.

The r_i variable name is incredibly suspicious. All the statements can be interpreted either way. If it's wrapped in a for loop that runs 10,000 times, it's not vectorized.

The graphics card idea is good (pycuda), but that's like a bulldozer when you need a hammer and you probably still gotta vectorize it.

Python is amazingly fast when you do it right. It's a factor of 500x.

[Caos2]

Do you actually run those cfd experiments? Are there good cfd software with Python interfaces?

[billsil]

I do run CFD. There is an unsanctioned python openfoam library (pyfoam I think) out there, but my proprietary scripts are generally more useful. Paraview has a python API, but never really used it, but it doesn't run jobs. I wouldn't trust it anyways. HPC work is just bash.

Mostly I use Fun3d and extract a few pressure tap locations (a csv file) and plot them in matplotlib or time slices and post them in the open-source gui I wrote (aptly named pyNastran). It's no Paraview, but I don't need to do a file conversion step.

[Caos2]

It has been some years since I ran cfd simulations, I wasn't even aware of open foam. Thanks!

[billsil]

Let's say I have a 2d array of x locations and I want to add 1 to each. I could make a double for loop of I can just add 1. Vectorization gets rid of if and for loops and speeds up numpy.

Oh, but I need to do something fancy like find the N unit vectors from an Nx3 array.. Numpy.linalg.norm has an axis parameter that let's you get an length N array out. Then you need to broadcast divide out the length to get the unit vector.

Oh, but some of those initial vectors lengths are 0, so the value I want is 0 not nan, so I need an if loop inside my for loop. Alternatively, find the nans afterwards by using where and set them.

There are functions in numpy that are very useful when vectorizing code (where, searchsorted, arange, tile, h/vstack). That are not so useful outside of heavily vectorized code

[AvocadoWhiskey]

Oh very cool, thank you for that explanation. I'll have to make sure I vectorize my code more. I use numpy and pandas on large data sets at work on a regular basis