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)
Improve this question
4
  • 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... Commented Nov 24, 2014 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)? Commented Nov 24, 2014 at 13:51
  • Done. It's really just a matrix on a vector field, and the results visualized by quivers arrows Commented Nov 24, 2014 at 14:03
  • I posted an answer on doing this kind of dot product here, may be helpful to you. Commented Nov 24, 2014 at 14:07

4 Answers 4

Reset to default
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)
Improve this answer
Sign up to request clarification or add additional context in comments.

2 Comments

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.
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()])
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)
Improve this answer

Comments

1

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.]])
Improve this answer

2 Comments

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 That's true, it was just because the OP was operating on ndarrays. 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)." docs.scipy.org/doc/numpy/reference/generated/…
0

If A is a matrix (2D array) and V is an array of vectors, then what you want is (A @ V[..., np.newaxis])[..., 0]. This works because [..., np.newaxis] turns each vector in V into a column vector (i.e. an n x 1 2D array).

If you're starting with an array of x values X and an array of y values Y, then you can turn them into an array of vectors [x, y] by setting V = np.column_stack([X, Y]).

Improve this answer

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.