3

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

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

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

This way, M is a representation of the function f over the cartesian product "v1 x v2", since M[i,j] = f(v1[i],v2[j]).

But this works because both sums and exponential work in parallel componentwise. My question is:

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

Please let me know if I should include more details!

rod
  • 141
  • 5
  • 1
    `np.exp(x1, x2)`? `numpy.exp(xx + yy)`? Which one? – Mad Physicist Oct 13 '22 at 14:42
  • Why are you using `meshgrid`, when you can use broadcasting? – Mad Physicist Oct 13 '22 at 14:44
  • Thanks, that was a typo, I even double checked twice! You mean the function numpy.broadcast? I didn't know it: I'm checking it! – rod Oct 14 '22 at 11:30
  • Ok, I read the "broadcasting" documentation page, but I wouldn't know how to exploit it – rod Oct 14 '22 at 12:26
  • I'm still not sure what you're trying to do. You have a NxN result from `x.T @ Q @ x`. Why do you need meshgrid? – Mad Physicist Oct 14 '22 at 15:31
  • To help clarify, make a couple of small arrays, say length 4 or 5, or a 2x4 array or so. Make a simple Q, like np.eye(2), and show what operations you want with the expected result. Right now, I'm having a hard time figuring out what your actual question is. – Mad Physicist Oct 14 '22 at 15:34
  • I am looking for a way to calculate `x.T @ Q @ x` on several inputs `x` at once. – rod Oct 14 '22 at 17:05
  • possible duplicate https://stackoverflow.com/questions/26849910/numpy-matrix-multiplication-broadcast – Him Oct 14 '22 at 17:16
  • @rod. Gotcha. Can you make that example then? I'd like to understand better what your inputs and expected outputs look like. If you can show me that, I suspect that I can make an answer for you. – Mad Physicist Oct 14 '22 at 19:26
  • @MadPhysicist I cannot make an example because I don't know how to do it! That is my question ^^ – rod Oct 18 '22 at 15:20
  • You absolutely can. Write a line that shows a sample input, and a hand calculated expected output. I understand that you want help vectorizing your operation, but showing the actual results you want for a small toy input will go a long way in making your prose understandable. If you can add a reference implementation using `for` loops, that would remove any remaining ambiguity. As it stands, I still don't have enough information to visualize how you imagine this happening. – Mad Physicist Oct 18 '22 at 15:37

1 Answers1

3
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 produces 2 arrays of a shape (len(y), len(x)), where first one is with x values along the second dimension. If we apply to this pair np.array then a 3d array of shape (2, len(y), len(x)) will be produced. With transpose we obtain an array, where an element indexed by [i,j,k] is x[i] if k==0 else y[j], where k is 0 or 1, i.e. first or second array from meshgrid.

With 'ijk,kn,ijn->ij' we tell einsum to return the sum written bellow for 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]].

Vitalizzare
  • 4,496
  • 7
  • 13
  • 32
  • Does `einsum` exploit nice matrix-multiplication algorithms here? I really love the `einsum` syntax, but IIRC it's not always the fastest `np` solution if one can use matrix ops. – Him Oct 14 '22 at 17:14
  • @Him I'm not sure if I understand your question. Actually, I don't know what is _nice matrix-multiplication algorithms_ and how to compare them with other solutions. Could you please explain what do you mean? Maybe you have a link to an article about this. – Vitalizzare Oct 14 '22 at 18:02
  • The way that you learn to multiply matrices "by hand" in school is not very efficient. Numerical libraries (e.g. `numpy`) all implement better algorithms for matrix multiplication. AFAIK, this doesn't extend to `einsum` because it is too "flexible" in some sense. https://en.wikipedia.org/wiki/Computational_complexity_of_matrix_multiplication#Simple_algorithms – Him Oct 14 '22 at 18:22
  • @Him The best alternative I can think of is `(a * np.inner(a, q)).sum(2)`. But on my machine it's much slower then `einsum`. What would you suggest as a replacement for einsum? – Vitalizzare Oct 14 '22 at 20:20
  • Thanks for your anwer Vitalizzare . I was going to ask what @Him asked. I.e.: is it efficient? – rod Oct 18 '22 at 15:22
  • 2
    np.einsum is generally more efficient than regular multiplication – Mad Physicist Oct 19 '22 at 07:20