Option to lambdify to use array instead of matrix

[asmeurer]

Issues like this seem to come up alot, where lambdify doesn't work like people expect because it converts Matrix to np.matrix instead of np.array. There should be some option to lambdify to make this easier to change, since the solution, lambdify(vars, expr, [{'ImmutableMatrix': np.array}, "numpy"]), is far from obvious.

Or maybe it can somehow be smarter in this case. I'm not sure how, though. As far as I can tell, it's just np.matrix that refuses to do things that work on a 2x2 array.

[moorepants]

Here's my reference issue on this subject: #2826

Yes, a very common use of lambdify it to transform SymPy matrices to numpy arrays to take advantage of vectorized functions. In fact, most everything in NumPy is based around the array not the matrix. Ideally we'd have SymPy arrays and SymPy matrices just like NumPy does. Then you create a SymPy array in Sympy if you want to eventually lambdify it with numpy to get an array.

[jcrist]

This has bugged me for a while now too. Would a default as np.array instead of np.matrix make sense? As @moorepants pointed out, most numpy operations focus on the array, so this as a default may make more sense.

On the other hand, sympy matrices define multiplication to be matrix multiplication. Having the original sympy matrix support A*B, while the lambdifed matrix supports A.dot(B) may be unexpected for users. Perhaps a kwarg boolean to define which behavior is wanted?

[moorepants]

If you read this potential PEP from the NumPy devs: numpy/numpy#4351 they declare that the matrix class in NumPy is essentially on its way out and is pretty much deprecated internally (maybe not for SciPy sparse). So it may be the right move, looking forward, to make this default to array.

But, I'm not sure SymPy even has the ability to create N-dimensional objects. If all we have are 2D arrays (matrices) in SymPy, then it does make sense for them to map that way.

As an intermediate solution, I am in favor of having a simpler kwarg to enable np.array to be used instead of np.matrix. It would keep things backwards compatible and not be too much of a burden for users. We'd need nice docs explaining this so it is quickly obvious when using looking at the docstring of lambdify.

[jcrist]

I added a kwarg "use_array". It defaults to False, which keeps the behavior of lambdify unchanged. Enabling it substitutes np.array for np.matrix. Should be a decent interim solution.

One thing that I wasn't sure how to handle was the case of a 1 col/ 1 row Sympy Matrix. In Sympy, matrices are always 2 dimensional. Some scipy and numpy functions only allow 1-dimensional arrays. For example, the scipy function odeint takes a function func which must return a 1-d object. Using the fix I just submitted, the workflow for this is the following:

x, y = symbols('x, y')
dx = lambdify((x, y), Matrix([y, (1-x**2)*y - x]), use_array=True)
func = lambda x, t: dx(*x).reshape(2)
x = odeint(func, array([0, 1]), linspace(0, 5, 1000))

The lambdify function could be made smarter to convert a 1 col/1 row Sympy matrix into a 1-d numpy array to remove the need for the reshape call, but I'm not sure if this is desirable in all situations.