Opened 3 years ago
Closed 21 months ago
#14403 closed defect (fixed)
Symbolic charpoly broken
Reported by: | nbruin | Owned by: | burcin |
---|---|---|---|
Priority: | major | Milestone: | sage-6.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 sage-devel.
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 3 years ago by eviatarbach
- Cc eviatarbach added
Changed 2 years ago by gagern
comment:2 Changed 2 years ago by gagern
comment:3 Changed 2 years ago by gagern
- Cc gagern added
comment:4 Changed 2 years ago by ppurka
- Cc ppurka added
comment:5 follow-up: ↓ 6 Changed 2 years ago by 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.
comment:6 in reply to: ↑ 5 ; follow-up: ↓ 23 Changed 2 years ago by ppurka
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 follow-ups: ↓ 8 ↓ 9 Changed 2 years ago by nbruin
Please look at the sage-devel 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 2 years ago by ppurka
Replying to nbruin:
Please look at the sage-devel 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 sage-devel 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 ; follow-up: ↓ 10 Changed 2 years ago by gagern
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 well-known 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 2 years ago by gagern
comment:10 in reply to: ↑ 9 ; follow-up: ↓ 11 Changed 2 years ago by nbruin
Replying to gagern:
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.
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 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.
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 an x. 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 ; follow-up: ↓ 13 Changed 2 years ago by gagern
Replying to nbruin:
The routine polynomial tries to recognize the thing as a polynomial in whatever variables (hence the problems with y^(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 of x, 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 call charpoly directly, but some other computation uses it.
- If you want to have an automatically chosen and non-conflicting 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 2 years ago by gagern
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 ; follow-up: ↓ 14 Changed 2 years ago by nbruin
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 of x, 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 call charpoly directly, but some other computation uses it.
Exactly. I agree.
- If you want to have an automatically chosen and non-conflicting 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 non-conflicting 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 ; follow-up: ↓ 15 Changed 2 years ago by gagern
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 non-conflicting 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 ; follow-up: ↓ 16 Changed 2 years ago by nbruin
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 low-cost, 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 though charpoly 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 2 years ago by gagern
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 2 years ago by burcin
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 2 years ago by gagern
Changed 2 years ago by gagern
comment:18 Changed 2 years ago by gagern
- Status changed from new to needs_review
comment:19 Changed 2 years ago by nbruin
testbot apply trac_14403_symbolic_charpoly_v3.patch trac_14403_expression_polynomial.patch
comment:20 Changed 2 years ago by jdemeyer
- Milestone changed from sage-5.11 to sage-5.12
comment:21 Changed 2 years ago by egourgoulhon
- Cc egourgoulhon added
comment:22 Changed 22 months ago by vbraun_spam
- Milestone changed from sage-6.1 to sage-6.2
comment:23 in reply to: ↑ 6 Changed 22 months ago by rws
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 22 months ago by rws
- Branch set to u/rws/ticket/14403
- Created changed from 04/02/13 22:15:26 to 04/02/13 22:15:26
- Modified changed from 02/18/14 03:02:47 to 02/18/14 03:02:47
comment:25 Changed 21 months ago by git
- 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 21 months ago by rws
- 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 21 months ago by vbraun
please fill in authors/reviewers.
comment:28 Changed 21 months ago by rws
- Reviewers set to Ralf Stephan
comment:29 Changed 21 months ago by vbraun
- 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 ready-to-use 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 letter x 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 the expression.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.