0

I'm trying to efficiently evaluate multiple instances of a python function that takes in 4 arguments, but I can't find how to map the function to several input values in a way that doesn't involve native python's slow for loops.

Here is an example of where I call the function just once (i.e. for one observation). 

from scipy.stats.mvn import mvnun

lower_bounds_one_obs = np.full(5, -1)
upper_bounds_one_obs = np.full(5,  1)

cov_mtx = np.eye(5)
means = np.zeros(5)

single_result = mvnun(lower_bounds_one_obs, upper_bounds_one_obs, means, cov_mtx)[0]

print(single_result)

For the sake of context, the mvnun function calculates the rectangular CDF of a multivariate normal distribution between lower and upper bounds. In this specific case, the lower and upper bounds have 5 elements in them.

What I want to do is to calculate mvnun for an array of lower and upper bounds, but always using the same means and cov_mtx.

In other words, I want to have two main input arrays, lower_bounds and lower_bounds, both with shape (1000 x 5). Then I would go through the process 1000 times, calculating mvnun for lower_bounds[0] and upper_bounds[0], then lower_bounds[1] and upper_bounds[1], and so on and so forth until lower_bounds[999] and upper_bounds[999].

In scipy's documentation, I found that I could vectorize/map this function, but only for the upper bounds, making the lower bound fixed. This can be done like so:

lower_fixed = np.full(5, -np.inf)
upper_bounds = np.random.rand(5 * 1000).reshape((1000, 5))
func1d = lambda upper_slice: mvnun(lower_fixed, upper_slice, means, cov_mtx)[0]
out = np.apply_along_axis(func1d, -1, upper_bounds)

This way, mvnun is evaluated multiple times with a fixed lower bound and iterating over all 1000 rows of upper_bound.

Does anyone know how to adjust the snippet of code above so that I can use an array of lower bounds instead of fixed lower bounds?

Thanks!

PS: This question is a bit of an update on this. There are also some interesting insights here too, but the solutions they present are pretty slow because they use regular python loops. More ideas here too.

Felipe D.
  • 1,157
  • 9
  • 19
  • 1) numpy.vectorize or things like lambda upper_slice np.apply_along_axis(...) doen't improve performance if that is what you wan't This things call Python function in a Python for loop. They are meant for convinience and not for performance. 2) Your Test-example is limited by the compiled mvdst.f, the wrapper is insignificant here. Solutions: use Intel-scipy (I encountered a 3x -speedup) and some multiprocessing, or wrap the fortran subroutine yourdelf (fast Fortran compiler eg. ifort needed) and write a implementation using cython (quite some work) – max9111 Jul 24 '18 at 08:05

1 Answers1

0

For reference only. Not sure if it's efficient enough. On my machine it's a little faster than np.fromiter combined with python map.

lower_bounds = np.full((1000,5), -np.inf)
upper_bounds = np.random.rand(1000, 5)

bounds = np.concatenate([lower_bounds, upper_bounds], axis=1)  
func1d = lambda bound: mvnun(bound[:5], bound[5:], means, cov_mtx)[0]
out = np.apply_along_axis(func1d, -1, bounds)
Sacry
  • 535
  • 4
  • 9