Opened 7 years ago

Last modified 5 years ago

#12101 needs_work defect

infinite recursion with exp on sparse matrix

Reported by: benjamin.peterson Owned by: jason, was
Priority: major Milestone: sage-6.4
Component: linear algebra Keywords:
Cc: dsm, rbeezer Merged in:
Authors: Karl-Dieter Crisman Reviewers: Burcin Erocal
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by kcrisman)

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.

Attachments (1)

trac_12101-sparse-matrix-exp.patch (1.1 KB) - added by kcrisman 7 years ago.

Download all attachments as: .zip

Change History (10)

comment:1 Changed 7 years ago by kcrisman

  • Summary changed from infinite recursion with exp on integer matrix to infinite recursion with exp on sparse matrix

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.

Changed 7 years ago by kcrisman

comment:2 Changed 7 years ago by kcrisman

  • Authors set to Karl-Dieter Crisman
  • Cc dsm rbeezer added
  • Description modified (diff)
  • Status changed from new to needs_review

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

comment:3 follow-up: Changed 6 years ago by burcin

  • Reviewers set to Burcin Erocal

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?

comment:4 in reply to: ↑ 3 Changed 6 years ago by 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.

comment:5 Changed 6 years ago by kcrisman

  • Status changed from needs_review to needs_work

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.

comment:6 Changed 6 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:7 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:8 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:9 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4
Note: See TracTickets for help on using tickets.