Improvements to lambdify for matrices

[moorepants]

I'd like to be able to create a sympy matrix of expressions, lambdify it, and then pass in and n x m numpy array of the m arguments to get back an array of n matrices.

In [50]: import sympy

In [51]: from sympy.abc import a, b, c, d

In [52]: expr = a**2 + sympy.cos(d) - b * a + (c**d / sympy.sin(a))

In [53]: mat_expr = sympy.Matrix([[expr, expr], [expr, expr]])

In [54]: import numpy as np

In [55]: an = np.random.random(5)

In [56]: bn = np.random.random(5)

In [57]: cn = np.random.random(5)

In [58]: dn = np.random.random(5)

In [63]: f = sympy.lambdify((a, b, c, d), mat_expr, modules='numpy')

In [64]: args = (an, bn, cn, dn)

In [65]: f(*args)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-65-528d283b8125> in <module>()
----> 1 f(*args)

/usr/local/lib/python2.7/dist-packages/numpy/__init__.pyc in <lambda>(a, b, c, d)

/usr/local/lib/python2.7/dist-packages/numpy/matrixlib/defmatrix.pyc in __new__(subtype, data, dtype, copy)
    270         shape = arr.shape
    271         if (ndim > 2):
--> 272             raise ValueError("matrix must be 2-dimensional")
    273         elif ndim == 0:
    274             shape = (1, 1)

ValueError: matrix must be 2-dimensional

Expected result would be a 5 x 2 x 2 numpy array.

You can do this with a for loop:

In [75]: for i in range(5):
   ....:     f(an[i], bn[i], cn[i], dn[i])

But it is much too slow for long expressions and large n. It would be nice if this translated to a vectorized version on numpy's side.

This may just be a matter of generate code that uses numpy.array instead of numpy.matrix:

In [77]: def compute(a, b, c, d):
    return np.mat([[a + b, c * c], [d * c, b + b]])
   ....: 

In [78]: compute(an, bn, cn, dn)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-78-ecbae16015ee> in <module>()
----> 1 compute(an, bn, cn, dn)

<ipython-input-77-dfe80aa828d7> in compute(a, b, c, d)
      1 def compute(a, b, c, d):
----> 2     return np.mat([[a + b, c * c], [d * c, b + b]])
      3 

/usr/local/lib/python2.7/dist-packages/numpy/matrixlib/defmatrix.pyc in asmatrix(data, dtype)
     94 
     95     """
---> 96     return matrix(data, dtype=dtype, copy=False)
     97 
     98 def matrix_power(M, n):

/usr/local/lib/python2.7/dist-packages/numpy/matrixlib/defmatrix.pyc in __new__(subtype, data, dtype, copy)
    270         shape = arr.shape
    271         if (ndim > 2):
--> 272             raise ValueError("matrix must be 2-dimensional")
    273         elif ndim == 0:
    274             shape = (1, 1)

ValueError: matrix must be 2-dimensional

In [79]: def compute(a, b, c, d):
    return np.array([[a + b, c * c], [d * c, b + b]])
   ....: 

In [80]: compute(an, bn, cn, dn)
Out[80]: 
array([[[ 1.22660376,  0.22579537,  0.50625436,  1.30750868,  0.75817236],
        [ 0.75518298,  0.5985355 ,  0.43101343,  0.32604234,  0.03012157]],

       [[ 0.58882765,  0.48415567,  0.49814087,  0.20791578,  0.07610386],
        [ 0.83998489,  0.25576972,  0.18090873,  1.77232675,  1.06712591]]])

In [81]: compute(an, bn, cn, dn).shape
Out[81]: (2, 2, 5)

[moorepants]

Shouldn't I be able to specify that SymPy MutableDenseMatrices become the numpy array function?

In [1]: from sympy import cos, sin, Matrix, lambdify

In [2]: from sympy.abc import a, b, c, d

In [3]: import numpy as np

In [4]: expr = a**2 + cos(d) - b * a + (c**d / sin(a))

In [5]: mat_expr = Matrix([[expr, expr], [expr, expr]])

In [6]: type(mat_expr)
Out[6]: sympy.matrices.dense.MutableDenseMatrix

In [7]: f = lambdify((a, b, c, d), mat_expr, modules=({'MutableDenseMatrix':np.array}, 'numpy'))

In [9]: an = np.random.random(5)

In [10]: bn = np.random.random(5)

In [11]: cn = np.random.random(5)

In [12]: dn = np.random.random(5)

In [13]: args = (an, bn, cn, dn)

In [14]: f(*args)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-14-528d283b8125> in <module>()
----> 1 f(*args)

/usr/local/lib/python2.7/dist-packages/numpy/__init__.pyc in <lambda>(a, b, c, d)

/usr/local/lib/python2.7/dist-packages/numpy/matrixlib/defmatrix.pyc in __new__(subtype, data, dtype, copy)
    270         shape = arr.shape
    271         if (ndim > 2):
--> 272             raise ValueError("matrix must be 2-dimensional")
    273         elif ndim == 0:
    274             shape = (1, 1)

ValueError: matrix must be 2-dimensional

[moorepants]

This worked! Cool. Not sure why I had to specify ImmutableMatrix though.

In [19]: f = lambdify((a, b, c, d), mat_expr, modules=({'ImmutableMatrix':np.array}, 'numpy'))

In [20]: f(*args)
Out[20]: 
array([[[ 1.00808799,  1.57085225,  0.45757548,  3.29441973,  2.71774198],
        [ 1.00808799,  1.57085225,  0.45757548,  3.29441973,  2.71774198]],

       [[ 1.00808799,  1.57085225,  0.45757548,  3.29441973,  2.71774198],
        [ 1.00808799,  1.57085225,  0.45757548,  3.29441973,  2.71774198]]])

[asmeurer]

I think it sympifies everything, converting all matrices to immutable matrices.

[moorepants]

I'll close this, as it is possible to do what I want with the current options.

[asmeurer]

I'm told that arrays are better in numpy than matrices, and that the only real reason to use matrices is for the syntactic sugar for multiplication. So maybe we should default to array.