6

I'm having a hard time getting into numpy. What I want in the end is a simple quiver plot of vectors that have been transformed by a matrix. I've read many times to just use arrays for matrices, fair enough. And I've got a meshgrid for x and y coordinates

X,Y = np.meshgrid( np.arange(0,10,2),np.arange(0,10,1) )
a = np.array([[1,0],[0,1.1]])

But even after googling and trying for over two hours, I can't get the resulting vectors from the matrix multiplication of a and each of those vectors. I know that quiver takes the component length as input, so the resulting vectors that go into the quiver function should be something like np.dot(a, [X[i,j], Y[i,j]]) - X[i,j] for the x-component, where i and j iterate over the ranges.

I can of course program that in a loop, but numpy has sooo many builtin tools to make these vectorized things handy that I'm sure there's a better way.

edit: Alright, here's the loop version.

import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(10,10))

n=10
X,Y = np.meshgrid( np.arange(-5,5),np.arange(-5,5) )
print("val test", X[5,3])
a = np.array([[0.5,0],[0,1.3]])
U = np.zeros((n,n))
V = np.zeros((n,n))
for i in range(10):
    for j in range(10):
        product = np.dot(a, [X[i,j], Y[i,j]]) #matrix with vector
        U[i,j] = product[0]-X[i,j]  # have to substract the position since quiver accepts magnitudes
        V[i,j] = product[1]-Y[i,j]

Q = plt.quiver( X,Y, U, V)
YXD
  • 31,741
  • 15
  • 75
  • 115
Basti
  • 2,228
  • 1
  • 19
  • 25
  • What do you expect the result to look like? `a.shape = (2,2), X.shape = (10,5), Y.shape = (10,5)`. I don't get your point... – jkalden Nov 24 '14 at 13:46
  • Can you show a simple version in a loop (to make it easier to work out what you are trying to do)? – Lee Nov 24 '14 at 13:51
  • Done. It's really just a matrix on a vector field, and the results visualized by quivers arrows – Basti Nov 24 '14 at 14:03
  • I posted an answer on doing this kind of dot product [here](http://stackoverflow.com/a/22081723/553404), may be helpful to you. – YXD Nov 24 '14 at 14:07

3 Answers3

6

You can either do matrix multiplication "manually" using NumPy broadcasting like this:

import numpy as np
import matplotlib.pyplot as plt

X,Y = np.meshgrid(np.arange(-5,5), np.arange(-5,5))
a = np.array([[0.5, 0], [0, 1.3]])

U = (a[0,0] - 1)*X + a[0,1]*Y
V = a[1,0]*X + (a[1,1] - 1)*Y

Q = plt.quiver(X, Y, U, V)

or if you want to use np.dot you have to flatten your X and Y arrays and combine them to appropriate shape as follows:

import numpy as np
import matplotlib.pyplot as plt

X,Y = np.meshgrid(np.arange(-5,5), np.arange(-5,5))
a = np.array([[0.5, 0], [0, 1.3]])

U,V = np.dot(a-np.eye(2), [X.ravel(), Y.ravel()])

Q = plt.quiver(X, Y, U, V)
Ondro
  • 997
  • 5
  • 8
  • Humm, I appreciate that. But that are both extremely limited cases, right? First is a hand-written matrix multiplication. That's pretty hacky and won't work for any higher dimensions etc. And the second is relying in the zeroes. – Basti Nov 24 '14 at 17:13
  • Both cases are for 2D. For 3D it should be something like `U,V,W = np.dot(a-np.eye(3), [X.ravel(), Y.ravel(), Z.ravel()])` – Ondro Nov 24 '14 at 17:22
5

As the docs says, for multidimensional data np.mul (or @) works in the following way:

For N dimensions it is a sum product over the last axis of a and the second-to-last of b:

dot(a, b)[i,j,k,m] = sum(a[i,j,:] * b[k,:,m])

This is not what we want here. However, there are simple alternatives that does not involve flattening-unflattening or manual matrix multiplication: np.tensordot and np.einsum.

The first is an example taken directly from the docs:

to do a matrix-matrix product with the left-most indices instead of rightmost, you can do np.einsum('ij...,jk...->ik...', a, b).

import numpy as np
import matplotlib.pyplot as plt

X,Y = np.meshgrid(np.arange(-5,5), np.arange(-5,5))
a = np.array([[0.5, 0], [0, 1.3]])

U, V = np.einsum('ij...,jk...->ik...', a - np.eye(2), np.array([X, Y]))
Q = plt.quiver(X, Y, U, V)

The second is application of simple np.tensordot. We just teach it to sum over second axis (colums) for the first argument and over first axis (rows) for the first argument.

import numpy as np
import matplotlib.pyplot as plt

X,Y = np.meshgrid(np.arange(-5,5), np.arange(-5,5))
a = np.array([[0.5, 0], [0, 1.3]])

U, V = np.tensordot(a - np.eye(2), np.array([X, Y]), axes=(1, 0))
plt.quiver(X, Y, U, V)
Community
  • 1
  • 1
Ilya V. Schurov
  • 7,687
  • 2
  • 40
  • 78
0

I have struggled with the same question and ended up using the numpy.matrix class. Consider the example below.

import numpy as np

>>> transformation_matrix = np.array([(1, 0, 0, 1),
...                                   (0, 1, 0, 0),
...                                   (0, 0, 1, 0),
...                                   (0, 0, 0, 1)])
>>> coordinates = np.array([(0,0,0),
...                         (1,0,0)])
>>> coordinates = np.hstack((coordinates, np.ones((len(coordinates), 1))))
>>> coordinates
array([[ 0.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  0.]])

In this case the numpy.matrix class helps. The following code gives the expected result by transposing the coordinates into column vectors and the designated matrix-multiplication overload of the numpy.matrix class.

>>> (np.asmatrix(transformation_matrix) * np.asmatrix(coordinates).T).T
matrix([[ 1.,  0.,  0.,  1.],
        [ 2.,  0.,  0.,  1.]])
user1556435
  • 966
  • 1
  • 10
  • 22
  • 1
    Why did you define them first as arrays, only to convert them to matrices later? You could've done `transformation_matrix = np.matrix([(1,0,0,1), (0,1,0,0), (0,0,1,0), (0,0,0,1)])` – R. Navega Jan 08 '19 at 05:07
  • 2
    @R.Navega That's true, it was just because the OP was operating on `ndarray`s. Note that `asmatrix()` does not make a copy of the data: "Unlike matrix, asmatrix does not make a copy if the input is already a matrix or an ndarray. Equivalent to matrix(data, copy=False)." https://docs.scipy.org/doc/numpy/reference/generated/numpy.asmatrix.html – user1556435 Jan 09 '19 at 09:47