4

Ultimately, I want to remove all explicit loops in the code below to take advantage of numpy vectorization and function calls in C instead of python.

Below is simplified for uses of numpy in python. I have the following quadratic function:

def quadratic_func(a,b,c,x):
    return a*x*x + b*x + c

I am trying to optimize choices of a,b,c given input data x and output data y of the same size (of course, this should be done by linear regression...but humor me). Say len(x)=100. Easy to vectorize with scalars a,b,c to get back a result of length 100.

Let's say that we know a,b,c should be inside of [-10,10] and I optimize by building a grid and picking the point with the min sum square error.

a=np.arange(-10.0, 10.01, 2.0)
nodes=np.array(np.meshgrid(a,a,a)).T.reshape(-1,3) #3-d cartesian product with array of nodes

For each of the 1331 nodes, I would like to calculate all 1331 of the length 100 return values. 

res=[]
x=np.random.uniform(-5.0,5.0, 100)
for node in nodes:
    res.append(quadratic_func(*node, x=x))

How can I take advantage of broadcasting so as to get my list of 1331 items each with 100 values that are the results of calling quadratic_func on x? Answer must use vectorization, broadcasting, etc to get the orders of magnitude speed improvements I am looking for. Also, the answer must use calls to quadratic_func - or more generally, my_func(*node, x=x).

In real life I am optimizing a non-linear function that is not even close to being convex and has many local minimums. It is a great functional form to use if I can get to the "right" local minimum - I already know how to do that, but would like to get there faster!

Divakar
  • 218,885
  • 19
  • 262
  • 358
FinanceGuyThatCantCode
  • 381
  • 3
  • 10

1 Answers1

1

One approach using a combination of broadcasting and np.einsum -

np.einsum('ij,jk->ik',nodes,x**np.array([2,1,0])[:,None])

Another one using matrix-multiplication with np.dot -

nodes.dot(x**np.array([2,1,0])[:,None])
Divakar
  • 218,885
  • 19
  • 262
  • 358
  • I see what you are doing, but looking for something more general than that - I don't want to have to know what the function is a priori. I want to be able to pass the objective function I am ultimately going to optimize parameters for into my optimizer and I can't assume that the functional form will be some combination of polynomials, sqrts, logs, etc. Can you do this using a call of quadratic_func? – FinanceGuyThatCantCode Oct 19 '16 at 21:31
  • @FinanceGuyThatCantCode Yeah that's the thing with NumPy vectorization, there's isn't a way that would take any generic func and speed it up. You can look into cython or numba if working with some specific functionality that could not be translated into NumPy funcs or ufuncs. – Divakar Oct 19 '16 at 21:34
  • I guess that is why I couldn't figure it out myself! Would be a nice feature though....and to think I could have asked the Continuum guys who were in my office today! – FinanceGuyThatCantCode Oct 19 '16 at 21:37