Danger of mixing numpy matrix and array

Question

The science/engineering application I'm working on has lots of linear algebra matrix multiplications, therefore I use Numpy matrices. However, there are many functions in python that interchangeably accept matrix or array types. Nice, no? Well, not really. Let me demonstrate the problem with an example:

from scipy.linalg import expm
from numpy import matrix

# Setup input variable as matrix
A = matrix([[ 0, -1.0,  0,  0],
            [ 0,  0,  0,  1.0],
            [ 0,  0,  0,  0],
            [ 0,  0,  1.0,  0]])

# Do some computation with that input
B = expm(A)

b1 = B[0:2, 2:4]
b2 = B[2:4, 2:4].T

# Compute and Print the desired output
print "The innocent but wrong answer:"
print b2 * b1

print "The answer I should get:"
print matrix(b2) * matrix(b1)

When run you get:

The innocent but wrong answer:
[[-0.16666667 -0.5       ]
 [ 0.          1.        ]]
The answer I should get, since I expected everything to still be matrices:
[[ 0.33333333  0.5       ]
 [ 0.5         1.        ]]

Any tips or advice on how to avoid this sort of a mix up? Its really messy to keep wrapping variables in matrix() calls to ensure they still are matrices. It seems there is no standard in this regard, and so it can lead to bugs that are hard to detect.

Answer 1

I tend to use array instead of matrix in numpy for a few reasons:

1. matrix is strictly 2D whereas you can have a numpy array of any dimension.
2. Aside from a few differences, array and matrix operations are pretty much interchangeable for a Matlab user.
3. If you use array consistently, then you would use numpy.dot() (or in Python 3.5 the new @ binary operator) for matrix multiplication. This will prevent the problem of not sure what * actually does in your code. And when you encounter a multiplication error, you can find the problem easier since you are certain of what kind of multiplication you are trying to perform.

So I would suggest you try to stick to numpy.array, but also keep in mind the differences between array and matrix.

Lastly, I found it a joy to work with numpy/scipy on bpython. The auto-prompt helps you to learn the properties of the function you are trying to use at a much faster pace than having to consult the numpy/scipy doc constantly.

Edit: The difference between array vs matrix is perhaps best answered here: 'array' or 'matrix'? Which should I use?

Answer 2

Mixing matrices and regular ndarrays can indeed be tricky and often not worth the hassle. I would second other posters and advise you to stick to arrays.

Nevertheless, in your particular example, the problem comes from expm. According to the doc, it takes a regular ndarray as argument and outputs a ndarray. If you want to transform your output back to matrix, you could use:

B = matrix(expm(A))

or

B = expm(A).view(matrix)

Now, B is a matrix, slices of B will be matrices themselves, and your multiplication will work as expected.

So, an advice would be to always check the type of the output of a function.