Opened 10 years ago
Last modified 18 months 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 )
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) ...
Attachments (1)
Change History (15)
comment:1 Changed 9 years ago by
- Summary changed from infinite recursion with exp on integer matrix to infinite recursion with exp on sparse matrix
Changed 9 years ago by
comment:2 Changed 9 years ago by
- 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: ↓ 4 Changed 9 years ago by
- 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:
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 9 years ago by
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 9 years ago by
- 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 8 years ago by
- Milestone changed from sage-5.11 to sage-5.12
comment:7 Changed 7 years ago by
- Milestone changed from sage-6.1 to sage-6.2
comment:8 Changed 7 years ago by
- Milestone changed from sage-6.2 to sage-6.3
comment:9 Changed 7 years ago by
- Milestone changed from sage-6.3 to sage-6.4
comment:10 Changed 18 months ago by
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.
comment:11 Changed 18 months ago by
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?
comment:12 follow-up: ↓ 13 Changed 18 months ago by
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.
comment:13 in reply to: ↑ 12 ; follow-up: ↓ 14 Changed 18 months ago by
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.
comment:14 in reply to: ↑ 13 Changed 18 months ago by
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.
This is actually a problem with sparse matrices (diagonal matrices are sparse). Here is an example.
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.
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.
I guess the answer would be to change this code to make sparse matrices dense?
Yup. Patch coming up.