Python – Define quadratic functions using numpy.meshgrid

Let’s consider a function of two variables f(x1, x2), where x1 spans vector v1 and x2 spans vector v2.

If f(x1, x2

) = np.exp(x1 + x2), we can represent this function as a matrix in Python by the command numpy.meshgrid as follows:

xx, yy = numpy.meshgrid(v1, v2)
M = numpy.exp(xx + yy)

Thus, M is the representation of the function f on the Cartesian product “v1 x v2” because M[i,j] = f(v1[i],v2[j]).

But this is possible because both the sum and the exponent are calculated in parallel. My question is:

If my variable is x = numpy.array([x1, x2]) and f is the quadratic function f(x) = x.T @ np .dot(Q, x), where Q is a 2×2 matrix, how can I do the same with the meshgrid function (i.e. calculate all values x of the function f on “v1″ v2” once)?

Please let me know if I should include more details!

def quad(x, y, q):
    """Return an array A of a shape (len(x), len(y)) with 
    values A[i,j] = [x[i],y[j]] @ q @ [x[i],y[j]]
    x, y: 1d arrays, 
    q: an array of shape (2,2)"""

from numpy import array, meshgrid, einsum
    a = array(meshgrid(x, y)).transpose()
    return einsum('ijk,kn,ijn->ij', a, q, a)

Notes

MeshGrid generates 2 arrays of shapes (len(y), len(x)), the first of which has x values along the second dimension. If we apply this pair to np.array, then a 3D array of shapes (2, len(y), len(x)) will be generated. Using Transpose we get an array where the elements indexed by [i,j,k] are x[i] if k==0 else y[j ], where k is 0 or 1, i.e. the first or second array from the meshgrid.

Using ‘ijk,kn,ijn->ij' we tell einsum to return the sum of each i, j:

sum(a[i,j,k]*q[k,n]*a[i,j,n] for k in range(2) for n in range(2))

Note that a[i,j]

== [x[i], y[j]].