Matrices and arrays in Python using np.array

Question

Let's say I have this calculation to do:

What is the best way to do it in Python? I come from a Matlab background and Matrices in Matlab are so easy to handle. I have tried with numpy but it gives me this error message:

ValueError: setting an array element with a sequence. 

Here is a snippet of my code:

for i in Xh:
    for u in Yh:    
        E= (np.array([(C,D),(E,F)]) * np.array([(i),(u)]) ) + np.array([Cx,Cy])

Please note that Xh and Yh are lists already calculated, this is why I am using the for loops (I know I can make it faster with list comprehensions).

Edit: This is a good solution that I got from this topic, and that works perfectly. For anyone who is facing the same problem, it could help.

a = np.matrix([[1,2],[3,4]])
v1 = np.matrix([[0],[1]])
v2 = np.matrix([[1],[1]])

res = np.dot(a,v1)+v2 # python version <= 3.5
res = a@v1+v2 # python version > 3.5 with the new dot operator @

Answer 1

If you want to have an experience close to matlab it is possible to use np.matrix instead of np.array. Your mathematical example would then be, e.g.

a = np.matrix([[1,2],[3,4]])
v1 = np.matrix([[0],[1]])
v2 = np.matrix([[1],[1]])

res = a*v1+v2

This is however not recommended by the numpy docs as pointed out in the comments. The way to do it with arrays is very similar. Here * is however the element-wise multiplication (.* in matlab). The matrix multiplication can then be done with:

a = np.matrix([[1,2],[3,4]])
v1 = np.matrix([[0],[1]])
v2 = np.matrix([[1],[1]])

res = np.dot(a,v1)+v2 # python version <= 3.5
res = a@v1+v2 # python version > 3.5 with the new dot operator @

Answer 2

I found two problems in your code:

1. You use round instead of squared brackets - you have to feed numpy's 2d arrays with a list of lists [[a0, a1],[a2, a3]]. This is where your ValueError comes from.
2. You use m1 * m2, which is the element-wise dot-product, but you should be using the matrix dot product np.dot(m1,m2) or its equivalent form m1.dot(m2).

Corrected:

for i in Xh:
    for u in Yh:
        E = np.dot(np.array([[C, D], [E, F]]), np.array([[i], [u]])) + np.array([[Cx], [Cy]])

Answer 3

You can make use of np.meshgrid to produce all combinations of Xh and Yh, and then np.einsum and broadcasting here to vectorize (read: speed up by a factor of 100x) the entire operation in one line:

import numpy as np

C, D, E, F = 1, 2, 3, 4
Cx, Cy = 1, 2
Xh = [1, 2, 3]
Yh = [4, 5, 6]

XhYh = np.array(np.meshgrid(Xh, Yh)).T.reshape(-1, 2)
M = np.array([[C, D], [E, F]])
V = np.array([Cx, Cy])

E = np.einsum('mk,fk->fm', M, XhYh) + V

Speed comparison for 100000 XhYh entries:

%timeit np.einsum('mk,fk->fm', M, XhYh) + V
# 209 µs ± 2.51 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit np.array([M.dot(np.array([x, y])) + V for x in Xh for y in Yh])
# 32.5 ms ± 608 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)