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.

[gh-pamaury]

Hi, this is still an issue with recent versions of sagemath (8.8). I think it would be a good idea to implement some proper error handling to tell the user what is going on. I don't really understand the code but printing our a message like

exp for sparse matrice is not implemented, please convert to a dense matrix using .dense_matrix()

would go a long way to make sage user friendly in this case.

[kcrisman]

Thanks, that is a good point. One could do that fairly easily - raise(NotImplementedError,"message") or the like. Would you be interested in trying to create such a branch?

[gh-mwageringel]

Well, this was fixed in #28935 (merged in 9.1.beta0) by converting to dense matrices, as I was unaware of this old ticket.

Perhaps, this ticket could be repurposed to implement exp for sparse matrices instead.

[kcrisman]

Well, this was fixed in #28935 (merged in 9.1.beta0) by converting to dense matrices, as I was unaware of this old ticket.

Perhaps, this ticket could be repurposed to implement exp for sparse matrices instead.

Yes, indeed! Please do. And/or to ask the user whether they want this ... unfortunately the comment about

converting to a dense matrix is something to be avoided

probably still applies.

[gh-mwageringel]

Replying to kcrisman:

Yes, indeed! Please do. And/or to ask the user whether they want this ... unfortunately the comment about

converting to a dense matrix is something to be avoided

probably still applies.

Personally, I am not that interested in this, as I think the change in #28935, although not ideal, is good enough for now. I think there are more important problems to solve regarding exp, for example the method currently only works for very small or structured matrices, so the conversion to a dense representation is not necessarily the bottleneck:

sage: matrix.random(ZZ, 4).exp()   # does not terminate in a reasonable amount of time
...
sage: matrix.random(ZZ, 5).exp()   # fails
...
TypeError: ECL says: Error executing code in Maxima: Unable to find the spectrum

Besides, the exponential of a sparse matrix is not necessarily sparse again, although sometimes it is:

sage: matrix.random(ZZ, 20, density=.25, sparse=True).change_ring(RDF).exp().density()
1

Computing the matrix exponential symbolically seems like a difficult problem in general.