Opened 7 years ago
Closed 6 years ago
#14403 closed defect (fixed)
Symbolic charpoly broken
Reported by:  nbruin  Owned by:  burcin 

Priority:  major  Milestone:  sage6.2 
Component:  symbolics  Keywords:  
Cc:  eviatarbach, gagern, ppurka, egourgoulhon  Merged in:  
Authors:  Martin von Gagern  Reviewers:  Ralf Stephan 
Report Upstream:  N/A  Work issues:  
Branch:  0434e05 (Commits)  Commit:  0434e05e991250f7eb0a19326be1a128a0baccce 
Dependencies:  Stopgaps: 
Description
see sagedevel. Not good:
sage: a=matrix([[x,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]) sage: charpoly(a) 0
Also not good
sage: var('y') sage: b=matrix([[y^(1/2),0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]) sage: charpoly(b) Error
Both can be solved by doing something more along the lines of
V=SR('crazy_varname') cp=self._maxima_(maxima).charpoly('crazy_varname')._sage_().expand() SR[var](list(cp.coefficient(V,i) for i in range(self.nrows()+1)))
in a.charpoly
.
Attachments (4)
Change History (33)
comment:1 Changed 7 years ago by
 Cc eviatarbach added
Changed 6 years ago by
comment:2 Changed 6 years ago by
comment:3 Changed 6 years ago by
 Cc gagern added
comment:4 Changed 6 years ago by
 Cc ppurka added
comment:5 followup: ↓ 6 Changed 6 years ago by
Thanks for taking on this ticket! It's a big problem.
There is in fact a way to generate unique variable names; just call SR.symbol()
. This will ensure that there are no collisions.
comment:6 in reply to: ↑ 5 ; followup: ↓ 23 Changed 6 years ago by
Replying to eviatarbach:
Thanks for taking on this ticket! It's a big problem.
There is in fact a way to generate unique variable names; just call
SR.symbol()
. This will ensure that there are no collisions.
That's true, but that generates ugly symbol names, which might be confusing to read in a long characteristic polynomial. I like the solution proposed in the patch a little better. For example, here is what I get from SR.symbol()
sage: SR.symbol() symbol160 sage: SR.symbol() symbol163 sage: SR.symbol() symbol166
Is there some general way of asking for the "next" symbol given some pattern? Essentially, emulating what a part of this patch does. Something which behaves like this:
sage: SR.next_variable('x') # suppose x, x0, x1 are already defined x2 sage: SR.next_variable('x2') x3 sage: SR.next_variable('x0y') # suppose that x0y is not yet defined. x0y
comment:7 followups: ↓ 8 ↓ 9 Changed 6 years ago by
Please look at the sagedevel thread referenced. Part of the problem here is the use of the routine "polynomial", which is dodgy for examples like:
matrix([[y^(1/2),0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])
You know that the expression returned by maxima's charpol is already presented as a polynomial in the (better unique!) variable used to do the expansion. The code proposed in the ticket is both more robust and faster.
comment:8 in reply to: ↑ 7 Changed 6 years ago by
Replying to nbruin:
Please look at the sagedevel thread referenced. Part of the problem here is the use of the routine "polynomial", which is dodgy for examples like:
matrix([[y^(1/2),0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])
You know that the expression returned by maxima's charpol is already presented as a polynomial in the (better unique!) variable used to do the expansion. The code proposed in the ticket is both more robust and faster.
Ok. I see. The main concern is this it seems (quoted from the sagedevel ML):
(You really don't want to let 'crazy_varname' leak out of charpoly, otherwise you'll have trouble if someone decides to make a matrix of characteristic polynomials and determine the charpoly of that. You also don't want to generate a fresh variable every time, since the management of variables in maxima/our interface would create a memory leak)
@gagern: can you try converting the code in the ticket description into a patch (along with a doctest)?
comment:9 in reply to: ↑ 7 ; followup: ↓ 10 Changed 6 years ago by
Replying to nbruin:
the use of the routine "polynomial" […] is dodgy […] The code proposed in the ticket is both more robust and faster.
I must confess I don't understand the difference. In my understanding, both the approach via coefficients
and the one via expression.polynomial
are intended to turn an expression representation of a polynomial into an element of the polynomial ring. I don't see why there is a difference in performance or robustness, but I can at least confirm a quite significant difference in performance. So I updated my patch accordingly, will attach that in a moment.
I did a mix of my old patch and the code mentioned in the ticket itself. In particular, I still use my own way to generate a unique variable. That way, chances are that different matrices using the same set of variables will use the same variable for the characteristic polynomial as well, which should avoid the memory leak scenario mentioned above although I don't fully understand that one either. Someone who does please confirm my assumption. Sounds like there were a central store of assigned variables, which I hadn't anticipated. Haven't looked at the SR sources yet.
In my setup, the variable name only depends on the matrix itself, not on any symbols used in other parts of the code, and more importantly not on the order in which characteristic polynomials are computed. I thought of trying a list of wellknown symbol names first, like ['x', 'y', 'z', 't', 'mu']
, but decided against it for now since it doesn't seem worth the effort.
There is no strong reason to use the same name in the maxima computation and for the symbolic ring. The original ticket used 'crazy_varname'
for one and var
for the other. On the other hand, I don't feel like relying on the fact that 'crazy_varname'
will be unlikely to occur in the matrix. I feel safer using a variable name which is certainly unique. And if I do the work to find such a vriable name, I might as well use the same for the computation in maxima and for the presentation as a polynomial. After all, always using x
as the variable name by default feels plain wrong if the matrix itself does contain an x
. That kind of overloading is bound to cause confusion.
My patch keeps the improvements for the expression.polynomial(…)
method. Even if charpoly
does not depend on it, being able to correctly convert a larger number of use cases seems worthwhile. If you'd prefer that as a separate commit, please let me know.
Changed 6 years ago by
comment:10 in reply to: ↑ 9 ; followup: ↓ 11 Changed 6 years ago by
Replying to gagern:
I must confess I don't understand the difference. In my understanding, both the approach via
coefficients
and the one viaexpression.polynomial
are intended to turn an expression representation of a polynomial into an element of the polynomial ring.
Just look up the source code of expression.polynomial
; you'll quickly see it's not a routine you want to call unless you really have to. The difference is that we already know with respect to which variable we want to collect the terms. The routine polynomial
tries to recognize the thing as a polynomial in whatever variables (hence the problems with y^(1/2)
)
Sounds like there were a central store of assigned variables, which I hadn't anticipated. Haven't looked at the SR sources yet.
Correct. This is the case in both maxima and in SR and necessarily permanently so in the translation layer between them (we cannot monitor the lifetime in both simultaneously)
There is no strong reason to use the same name in the maxima computation and for the symbolic ring. The original ticket used
'crazy_varname'
for one andvar
for the other. On the other hand, I don't feel like relying on the fact that'crazy_varname'
will be unlikely to occur in the matrix.
How about do_not_use_this_name_in_a_matrix_youll_compute_a_charpol_of
then? It saves you digging through all the variables occurring in the matrix unnecessarily.
I feel safer using a variable name which is certainly unique. And if I do the work to find such a vriable name, I might as well use the same for the computation in maxima and for the presentation as a polynomial. After all, always using
x
as the variable name by default feels plain wrong if the matrix itself does contain anx
. That kind of overloading is bound to cause confusion.
For me, explicit is better than implicit: If I'm going to do something with a charpol, I need to know in what variable it has been expressed, so the system choosing a relatively unpredictable name is no good.
For one thing, one might expect that
charpol(A)*charpol(B) == charpol(block_diagonal(A,B))
but that's not necessarily the case if sage varies the default variable name used.
comment:11 in reply to: ↑ 10 ; followup: ↓ 13 Changed 6 years ago by
Replying to nbruin:
The routine
polynomial
tries to recognize the thing as a polynomial in whatever variables (hence the problems withy^(1/2)
)
Not with my modification, but the speed difference remains. The reason why I'd have considered polynomial
to be potentially faster than the approach based on coefficient
is that the former has to traverse the expression only once, whereas the latter has to do one pass per power.
How about
do_not_use_this_name_in_a_matrix_youll_compute_a_charpol_of
then? It saves you digging through all the variables occurring in the matrix unnecessarily.
Possible for the computation, but we very certainly don't want this variable occuring in the result. Since I still maintain that we have to find a unique and readable name for the generator of the polynomial ring for the result, would there be any win in always using this name for the computation? Would this reduce memory leaks? I guess it might, so I'll include this in my next patch, but I'll wait for some answers to questions below first.
For me, explicit is better than implicit: If I'm going to do something with a charpol, I need to know in what variable it has been expressed, so the system choosing a relatively unpredictable name is no good.
I disagree. The explicit approach is always possible, i.e. you can always specify the variable name in the invocation. But if you only want to have a look at the polynomial, it doesn't matter that much what the variable is, as long as you can recognize it and can distinguish it from other variables. And if you want to do some subsequent computation, chances are that you don't care for the name of the variable at all, but only for the list of coefficients.
For one thing, one might expect that
charpol(A)*charpol(B) == charpol(block_diagonal(A,B))but that's not necessarily the case if sage varies the default variable name used.
I don't think I'd share that expectation. What I'd expect is
A.charpoly('t')*B.charpoly('t') == block_diagonal(A, B).charpoly('t')
i.e. if you explicitely name the variable, then the polynomials should agree. If you don't name the variable, and the default name of x
appears in one of the matrices but not the other, then raising an error instead would be acceptable to me, like what my patch would cause:
TypeError: unsupported operand parent(s) for '*': 'Univariate Polynomial Ring in x1 over Symbolic Ring' and 'Univariate Polynomial Ring in x over Symbolic Ring'
What exactly is it you propose for the visible result of this method, when called without a variable name?
 If you want a default of
x
in all cases, then you accept strange cases where the printed result will contain two flavors ofx
, one the generator of the polynomial and the other a variable from the symbolic ring. I consider this highly confusing. And it might cause even more trouble later on, particularly if this polynomial would get passed to some other system in text form, where the distinction would get lost.  If you want the
do_not_use…
name to leak into the result, reading the polynomial directly would be a real pain.  If you want to raise an exception whenever
x
occurs in the matrix, I'd consider this unneccessarily annoying. Particularly for cases where the user didn't callcharpoly
directly, but some other computation uses it.  If you want to have an automatically chosen and nonconflicting name, then this would be something along the lines I proposed, so any important difference in what you have in mind to what I did should be made more explicit.
comment:12 Changed 6 years ago by
I just noticed that charpoly
the function unconditionally passes its variable name to the method. Changed that, will be included in my next patch, once we have reached consensus on what the variable name in the default case should be. Of course if we were to agree that it should always be x
then that part would be superfluous, but I don't see me agreeing to that unless I see a really strong argument for it.
comment:13 in reply to: ↑ 11 ; followup: ↓ 14 Changed 6 years ago by
Replying to gagern:
What exactly is it you propose for the visible result of this method, when called without a variable name?
Indeed, x
, because that is already the default in sage anyway.
sage: P.<x>=QQ[] sage: M=matrix(2,2,[x,1,1,x]) sage: M.charpoly() x^2  2*x*x + x^2  1
 If you want a default of
x
in all cases, then you accept strange cases where the printed result will contain two flavors ofx
, one the generator of the polynomial and the other a variable from the symbolic ring. I consider this highly confusing. And it might cause even more trouble later on, particularly if this polynomial would get passed to some other system in text form, where the distinction would get lost.
You'd have to come up with a naming solution that works across sage structures, which sounds pretty hopeless to me.
 If you want the
do_not_use…
name to leak into the result, reading the polynomial directly would be a real pain.
It should definitely not leak! Then it can find its way back into a matrix and cause real havoc.
 If you want to raise an exception whenever
x
occurs in the matrix, I'd consider this unneccessarily annoying. Particularly for cases where the user didn't callcharpoly
directly, but some other computation uses it.
Exactly. I agree.
 If you want to have an automatically chosen and nonconflicting name, then this would be something along the lines I proposed, so any important difference in what you have in mind to what I did should be made more explicit.
I think choosing a nonconflicting name is too difficult and confusing (since the problem is not confined to just SR), so I'd accept a default that likely clashes as the best compromise. Conversion to other systems is indeed dangerous: don't do that :).
comment:14 in reply to: ↑ 13 ; followup: ↓ 15 Changed 6 years ago by
Replying to nbruin:
Indeed,
x
, because that is already the default in sage anyway.
I don't know sage development structures and politics. Is this a personal opinion, or does this come with enough authority to indicate this as the direction to go?
sage: P.<x>=QQ[] sage: M=matrix(2,2,[x,1,1,x]) sage: M.charpoly() x^2  2*x*x + x^2  1
Personally, I'd consider this a bug.
You'd have to come up with a naming solution that works across sage structures, which sounds pretty hopeless to me.
How many rings are there? How many of them have elements which you can name? While I agree that this might require considerable effort, I do have hope and would consider such an endeavour worthwhile.
I think choosing a nonconflicting name is too difficult and confusing
It seems our individual perceptions of how confusing the various alternatives are differ quite a lot… Can you point out an example where this confusion will actually cause harm?
Conversion to other systems is indeed dangerous: don't do that :).
This sounds both crucial and problematic to me. There is nothing to prevent users from doing this, and checking every possible conversion input for possible clashes sounds a lot more problematic than detecting symbol names in all structures which can carry them. So users will do this, and will get wrong results without a warning.
Besides, I was under the impression that the whole idea behind sage is that users don't have to think about what backend does the actual work, since it's all glued together nicely. Your statement sounds like a conscious violation of that ideal.
One thing which I can do is this: adapt my patch to do only the fix you quoted plus a few doctests. Then I would file another ticket for the fixes to the expression.polynomial
even though charpoly
no longer needs them, and yet another one for the broader question of choosing a suitable variable name for charpoly
(and similar situations if you can think of some). There we could continue to discuss the pros and cons of the two basic approaches: fixed but clashing vs. unpredictable unique. That way, the discussion about the correct approach to that could continue but would not block the required fix here. And if it were addressed eventually, it could be addressed in a way suitable to all affected structures. How does that sound?
comment:15 in reply to: ↑ 14 ; followup: ↓ 16 Changed 6 years ago by
Replying to gagern:
I don't know sage development structures and politics. Is this a personal opinion, or does this come with enough authority to indicate this as the direction to go?
As far as I'm concerned it's common sense. I won't claim authority for my own opinion concerning sage, though.
sage: P.<x>=QQ[] sage: M=matrix(2,2,[x,1,1,x]) sage: M.charpoly() x^2  2*x*x + x^2  1Personally, I'd consider this a bug.
I agree it's nasty, but it's there. The problem is that you can have arbitrarily nastily nested structures: QQ['s']['x']['u']['y']
. People could even name their number field generators or finite field generators x
.
How many rings are there? How many of them have elements which you can name? While I agree that this might require considerable effort, I do have hope and would consider such an endeavour worthwhile.
It has to be lowcost, though. You're competing with choosing x
just blindly.
One thing which I can do is this: adapt my patch to do only the fix you quoted plus a few doctests. Then I would file another ticket for the fixes to the
expression.polynomial
even thoughcharpoly
no longer needs them,
You might as well fix polynomial.
and yet another one for the broader question of choosing a suitable variable name for
charpoly
(and similar situations if you can think of some).
Yes, because I think it's not clear yet how to solve this properly and it's something that is a wider problem than just SR.
comment:16 in reply to: ↑ 15 Changed 6 years ago by
Replying to nbruin:
Replying to gagern:
and yet another one for the broader question of choosing a suitable variable name for
charpoly
(and similar situations if you can think of some).Yes, because I think it's not clear yet how to solve this properly and it's something that is a wider problem than just SR.
Filed #14972 for this. Will attach a patch for the rest shortly.
comment:17 Changed 6 years ago by
FWIW, I agree with Nils that it is""common sense" to use a standard variable name. There is no need to complicate the implementation by searching for variables.
If we decide to use a new variable name, the procedure to do this should be constant time. Enumerating variables in a matrix/ring and finding a new name is too complicated. Pick a prefix and increase a counter whenever you need a new variable. On the other hand, it would also be confusing to get different variables when you call charpoly()
twice.
Changed 6 years ago by
Changed 6 years ago by
comment:18 Changed 6 years ago by
 Status changed from new to needs_review
comment:19 Changed 6 years ago by
testbot apply trac_14403_symbolic_charpoly_v3.patch trac_14403_expression_polynomial.patch
comment:20 Changed 6 years ago by
 Milestone changed from sage5.11 to sage5.12
comment:21 Changed 6 years ago by
 Cc egourgoulhon added
comment:22 Changed 6 years ago by
 Milestone changed from sage6.1 to sage6.2
comment:23 in reply to: ↑ 6 Changed 6 years ago by
Replying to ppurka:
Is there some general way of asking for the "next" symbol given some pattern? Essentially, emulating what a part of this patch does. Something which behaves like this:
sage: SR.next_variable('x') # suppose x, x0, x1 are already defined x2 sage: SR.next_variable('x2') x3 sage: SR.next_variable('x0y') # suppose that x0y is not yet defined. x0y
I have opened #15831 for this. Also, below is the patch, adapted to base it on 6.0.
comment:24 Changed 6 years ago by
 Branch set to u/rws/ticket/14403
 Created changed from 04/03/13 05:15:26 to 04/03/13 05:15:26
 Modified changed from 02/18/14 11:02:47 to 02/18/14 11:02:47
comment:25 Changed 6 years ago by
 Commit set to 0434e05e991250f7eb0a19326be1a128a0baccce
Branch pushed to git repo; I updated commit sha1. New commits:
0434e05  Merge branch 'develop' into ticket/14403

comment:26 Changed 6 years ago by
 Status changed from needs_review to positive_review
The code looks good, matrix/ symbolics/ and misc/ long tests pass, sage docbuild works. So I think I'll set positive.
comment:27 Changed 6 years ago by
please fill in authors/reviewers.
comment:28 Changed 6 years ago by
 Reviewers set to Ralf Stephan
comment:29 Changed 6 years ago by
 Branch changed from u/rws/ticket/14403 to 0434e05e991250f7eb0a19326be1a128a0baccce
 Resolution set to fixed
 Status changed from positive_review to closed
OK, here is an attempt to address this.
I modified the
charpoly
method to always use a unique variable name. I found no readytouse utility function to achieve this, but I expect there might be, so if you know one please let me know and I'll update my code. So far I simply append numbers to the letterx
until the resulting string is not among the list of variable names used inside the matrix. The resulting variable name will always be used for the cached polynomial, and will by default be used for the returned polynomial.If the caller overrides the variable name, there is no check to ensure that this name does not conflict with the names used in the matrix. I'm unsure myself whether I'd consider this behavior a bug or a feature, i.e. what intended behavior for that case should be. If you want I can add a piece of code to perform that check before the
change_variable_name
call, and throw an error on conflicting variable names.I've kept the
expression.polynomial(None, ring=SR[var])
approach, and not used theexpression.coefficient(var, i)
version suggested in the ticket report. Instead, I adjusted the conversion of expressions to polynomials to only treat those variable powers as possible monomials where the variable itself corresponds to one of the variables from the target ring. All other expressions are treated as coefficients, even if they are powers of a symbolic variable. This avoids the conversion to integers for these cases, and hence the error mentioned in the ticket. This fix might be useful to other applications as well.