Numpy: How to vectorize parameters of a functional form of a function applied to a data set
<p>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.</p>

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

<pre><code>def quadratic_func(a,b,c,x):
    return a*x*x + b*x + c
</code></pre>

<p>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.</p>

<p>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.</p>

<pre><code>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
</code></pre>

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

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

<p>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).</p>

<p>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!</p>

<p>One approach using a combination of <a href="https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html" rel="nofollow"><code>broadcasting</code></a> and <a href="https://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html" rel="nofollow"><code>np.einsum</code></a> -</p>

<pre><code>np.einsum('ij,jk-&gt;ik',nodes,x**np.array([2,1,0])[:,None])
</code></pre>

<p>Another one using matrix-multiplication with <a href="https://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html" rel="nofollow"><code>np.dot</code></a> -</p>

<pre><code>nodes.dot(x**np.array([2,1,0])[:,None])
</code></pre>

python numpy vectorization numpy-broadcasting