CODEWITHSUNDEEP
numpyscipy
Basti
Basti

Reputation: 2448

Numpy: Multiply a matrix with an array of vectors

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)

Upvotes: 6

Views: 10395

Answers (4)

heptagon
heptagon

Reputation: 31

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]).

Upvotes: 0

Ilya V. Schurov
Ilya V. Schurov

Reputation: 8097

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)

Upvotes: 5

user1556435
user1556435

Reputation: 1056

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.]])

Upvotes: 1

Ondro
Ondro

Reputation: 1027

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)

Upvotes: 6

Related Questions