infinite recursion with exp on sparse matrix

[benjamin.peterson]

sage: exp(diagonal_matrix([1, 2, 3]))
  File "/sagenb/sage_install/sage-4.7.2/local/lib/python2.6/site-packages/sage/functions/log.py", line 126, in __call__
    dont_call_method_on_arg=dont_call_method_on_arg)
  File "function.pyx", line 715, in sage.symbolic.function.GinacFunction.__call__ (sage/symbolic/function.cpp:6666)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)
  File "matrix2.pyx", line 9933, in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:46796)

...

---

Apply trac_12101-sparse-matrix-exp.patch .

[kcrisman]

This is actually a problem with sparse matrices (diagonal matrices are sparse). Here is an example.

sage: D = matrix(SR,[1],sparse=True)
sage: type(D)
<type 'sage.matrix.matrix_generic_sparse.Matrix_generic_sparse'>
sage: D.exp()
ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (7994, 0))
<snip recursion>
/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/matrix/matrix2.so in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:52642)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/matrix/matrix2.so in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:52642)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/matrix/matrix2.so in sage.matrix.matrix2.Matrix.exp (sage/matrix/matrix2.c:52635)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/matrix/matrix_sparse.so in sage.matrix.matrix_sparse.Matrix_sparse.change_ring (sage/matrix/matrix_sparse.c:2220)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/matrix/matrix_generic_sparse.so in sage.matrix.matrix_generic_sparse.Matrix_generic_sparse.__copy__ (sage/matrix/matrix_generic_sparse.c:4362)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/matrix/matrix_generic_sparse.so in sage.matrix.matrix_generic_sparse.Matrix_generic_sparse.__init__ (sage/matrix/matrix_generic_sparse.c:2999)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/structure/element.so in sage.structure.element.Element.is_zero (sage/structure/element.c:6297)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.__nonzero__ (sage/symbolic/expression.cpp:10268)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.__nonzero__ (sage/symbolic/expression.cpp:10116)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.test_relation (sage/symbolic/expression.cpp:10932)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/rings/complex_interval_field.pyc in __call__(self, x, im)
    290 
    291             try:
--> 292                 return x._complex_mpfi_( self )
    293             except AttributeError:
    294                 pass

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._complex_mpfi_ (sage/symbolic/expression.cpp:6142)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._eval_self (sage/symbolic/expression.cpp:5539)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/symbolic/pynac.so in sage.symbolic.pynac.py_float (sage/symbolic/pynac.cpp:9642)()

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/rings/complex_interval_field.pyc in __call__(self, x, im)
    297             except AttributeError:
    298                 pass
--> 299         return complex_interval.ComplexIntervalFieldElement(self, x, im)
    300 
    301     def _coerce_impl(self, x):

/Users/.../sage-5.1.beta5/local/lib/python2.7/site-packages/sage/rings/complex_interval.so in sage.rings.complex_interval.ComplexIntervalFieldElement.__init__ (sage/rings/complex_interval.c:3364)()

RuntimeError: maximum recursion depth exceeded

If I make the zero matrix instead, it cuts off at the matrix init instead, not getting into the complex and pynac stuff.

I don't know why this is the particular thing that returns. I do know why we have an infinite recursion.

        from sage.symbolic.ring import SR
        return self.change_ring(SR).exp()

is the entire code in matrix2.pyx for the exponential method of a generic matrix. And the doctests only test dense matrices, whose coercion to the symbolic ring have nice exp methods. But sparse ones don't go anywhere.

sage: type(D)
<type 'sage.matrix.matrix_generic_sparse.Matrix_generic_sparse'>
sage: C = D.change_ring(SR)
sage: type(C)
<type 'sage.matrix.matrix_generic_sparse.Matrix_generic_sparse'>

I guess the answer would be to change this code to make sparse matrices dense?

sage: D.change_ring(SR).dense_matrix().exp()
[e]

Yup. Patch coming up.

[kcrisman]

Okay, this at least covers the basic case, has a couple tests, and checks that D is not changed by the operation in tests.

Cc:ing a couple reviewers who might be interested in matrices. I'd ask for any reviewer to try to find more obscure sparse matrices that might still not work; might as well take care of them at this time too. But most should have an exp method, I think, when made dense.

patchbot: Apply trac_12101-sparse-matrix-exp.patch

[burcin]

I don't think silently converting the sparse input matrix to a dense one is a good idea. We should define an exp() method for sparse symbolic matrices to avoid this infinite recursion.

Here is the code for the exp() method of Matrix_symbolic_dense:

    def exp(self):
        if not self.is_square():
            raise ValueError, "exp only defined on square matrices"
        if self.nrows() == 0:
            return self
        # Maxima's matrixexp function chokes on floating point numbers
        # so we automatically convert floats to rationals by passing
        # keepfloat: false
        m = self._maxima_(maxima)
        z = maxima('matrixexp(%s), keepfloat: false'%m.name())
        if self.nrows() == 1:
            # We do the following, because Maxima stupidly exp's 1x1
            # matrices into non-matrices!
            z = maxima('matrix([%s])'%z.name())

        return z._sage_()

It would be great if we could avoid calling maxima for this. How hard would it be to implement what maxima does natively in Sage? Here is the code for the matrixexp maxima function:

 http://maxima.git.sourceforge.net/git/gitweb.cgi?p=maxima/maxima;a=blob;f=share/linearalgebra/matrixexp.lisp;hb=HEAD

Another option is to find a way to convert a sparse matrix to Maxima and still use its matrixexp() implementation. Does Maxima have a sparse matrix constructor?

[tscrim]

Replying to burcin:

I don't think silently converting the sparse input matrix to a dense one is a good idea. We should define an exp() method for sparse symbolic matrices to avoid this infinite recursion.

My 2 cents; converting to a dense matrix is something to be avoided because (in general) it requires substantial memory allocation (ex. take a 1000x1000 matrix with 2 (non-zero) entries). Thus I would rather see the exp() implemented for sparse (symbolic) matrices and return a sparse matrix.

[kcrisman]

Of course, we're not actually converting to a dense matrix per se, we're just using a dense version of the matrix to do this (I hope). But if there is a way to keep things sparse, absolutely, that makes great sense.