Opened 5 years ago

Last modified 4 years ago

#16222 needs_info enhancement

Faster exactification using numeric minpoly

Reported by: gagern Owned by:
Priority: major Milestone: sage-6.4
Component: number fields Keywords:
Cc: Merged in:
Authors: Martin von Gagern Reviewers: Vincent Delecroix
Report Upstream: N/A Work issues:
Branch: u/vdelecroix/16222 (Commits) Commit: 8773468739544ed9f12b836c7f28a06be37ffd94
Dependencies: #15605 Stopgaps:

Description

This is a spinoff from comment:13:ticket:14239. There I noticed that for a large symbolic expressions b, the call b.minpoly() was a lot (by several orders of magnitude) faster than the call QQbar(b).minpoly(). We should try to make this speed gain available to those algebraic numbers which were constructed from symbolic expressions.

I now know that the speed gain is almost certainly due to the numeric algorithm in calculus.minpoly. So what we should in my opinion do is use that numeric algorithm to obtain a minimal polynomial of the whole expression, and if that succeeds, base subsequent exactifications on that polynomial instead of the nested tree of algebraic number descriptors.

Change History (26)

comment:1 Changed 5 years ago by gagern

  • Branch set to u/gagern/ticket/16222
  • Created changed from 04/23/14 19:33:11 to 04/23/14 19:33:11
  • Modified changed from 04/23/14 19:33:11 to 04/23/14 19:33:11

comment:2 Changed 5 years ago by gagern

  • Commit set to 491dcc4884823cab627f318ea63f0e2504f9f7c2

These modifications are not ready for inclusion yet, they need some discussion first.

  1. How to get the symbolic expression from expression_conversion into AlgebraicNumber_base.

Right now, I use direct access to an underscore-prefixed variable called _symbolic_expression. I'm not sure whether that is acceptable. If cross-module interfaces should avoid using underscores, we'll have to look for alternatives. But I rely on the fact that the expression really describes the number in question, so this certainly isn't a variable the user should be allowed to tweak. So I don't know a better way yet.

  1. Default state of that variable.

Right now, I have a class variable which defaults to None and instance variables which will override this when required. Alternatives would be setting an instance variable in some (as of now non-existant) constructor, or avoding a default value altogether and using hasattr instead. I simply wrote what I'd do in my own code, but if there is a sage convention to write things differently, I'll gladly comply.

  1. When to compute the minpoly.

I now do the computation on demand, hoping that expressions cannot be legally modified after their creation. If they can be modified, we might want to do the minpoly computation when we do the conversion from symbolic to algebraic. Likewise if you have reasons to believe that minpoly computation should not make that conversion noticably slower than it already is.

  1. Test case dependency.

I just noticed that my current doctest depends on ticket:14239, but that should be easy to fix. If you have a favorite example, let me know, but otherwise I'll come up with something.

  1. The _more_precision loop.

I have more trouble testing that _more_precision codepath inside _exact_symbolic. Can anyone come up with a situation where identifying the correct root would need more than default precision but where numeric minpoly computation would still succeed? If not, what shall we do? Shall we deactivate that whole loop, and simply return without taking action if we find more than one root which could match?

  1. How to obtain the roots.

Is p.roots(self.parent(), False) the right thing to do, or is this too high-level, should I use something more basic here?


New commits:

bb01e0bDrive-by fix to minpoly documentation.
491dcc4Faster exactification using numeric minpoly.

comment:3 Changed 5 years ago by git

  • Commit changed from 491dcc4884823cab627f318ea63f0e2504f9f7c2 to 8597313c855ce7a413cfcabead498059d2cfcbf6

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

1b703c0Drive-by fix to minpoly documentation.
8597313Faster exactification using numeric minpoly.

comment:4 Changed 5 years ago by gagern

  • Status changed from new to needs_review

Rebased to current develop, and changed test case to avoid dependency on #14239. I'm now requesting review since that's what I need even if I consider it likely that the reviewer will not agree with all the choices I made as outlined above.

comment:5 Changed 5 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4

comment:6 follow-up: Changed 5 years ago by vdelecroix

Hello,

Did you try using the magic algdep from pari? It perform (very quickly) some LLL to find the polynomial with smallest coefficient for which x is almost a root. With your example

sage: a = (443/96*I*sqrt(443)*sqrt(3) + 833939/1728)^(1/3)
sage: b = sqrt(144*a + 9205/a + 1176)/12
sage: c = QQbar(b)
sage: gp.algdep(c.n(prec=150).real(), 6)
16*x^6 - 392*x^4 + 133*x^2 + 900

The advantage is that it would potentially speed up non symbolic cases.

Vincent

comment:7 in reply to: ↑ 6 ; follow-up: Changed 5 years ago by gagern

Replying to vdelecroix:

Did you try using the magic algdep from pari?

My code builds on that, since expression.minpoly calls calculus.minpoly which in turn calls algdep. But it does come with a verification step, ensuring that the found polynomial is in fact the one we need. I'm not sure I'd trust the output from algdep directly, particularly not unless I had some very good ideas of what degree to expect.

The advantage is that it would potentially speed up non symbolic cases.

Well, one alternative would be to directly feed the descriptor DAG into calculus.minpoly. Computing a numerical_approx from these would be easy enough. The hard part would be the symbolic verification step, g(ex).simplify_trig().canonicalize_radical() == 0. I guess that's done in Maxima. Do you have an idea how we might either implement this for the non-symbolic case, or find a rigorous argument why we don't need such a proof?

Last edited 5 years ago by gagern (previous) (diff)

comment:8 in reply to: ↑ 7 ; follow-up: Changed 5 years ago by vdelecroix

Replying to gagern:

Replying to vdelecroix:

Did you try using the magic algdep from pari?

My code builds on that, since expression.minpoly calls calculus.minpoly which in turn calls algdep. But it does come with a verification step, ensuring that the found polynomial is in fact the one we need. I'm not sure I'd trust the output from algdep directly, particularly not unless I had some very good ideas of what degree to expect.

Of course we can not believe the output of algdep. But at least we can bound the degree. And once a polynomial is given we can check irreducibility. So I guess we can produce a candidate quite reliably.

The advantage is that it would potentially speed up non symbolic cases.

Well, one alternative would be to directly feed the descriptor DAG into calculus.minpoly. Computing a numerical_approx from these would be easy enough. The hard part would be the symbolic verification step, g(ex).simplify_trig().canonicalize_radical() == 0. I guess that's done in Maxima. Do you have an idea how we might either implement this for the non-symbolic case, or find a rigorous argument why we don't need such a proof?

Right, checking that it is actually an annihilator polynomial is harder... I do not see anything different from calling __nonzero__ which indeed calls .exactify().

comment:9 in reply to: ↑ 8 Changed 5 years ago by gagern

Replying to vdelecroix:

So I guess we can produce a candidate quite reliably.

I'm unsure about the relation between “candidate” and “reliable”. The former implies something that needs to be verified, the latter like I could rely on it without verification. But I guess what you mean is that we could construct a candidate which would very likely be the correct one. But “very likely” isn't enough in my opinion.

Right, checking that it is actually an annihilator polynomial is harder... I do not see anything different from calling __nonzero__ which indeed calls .exactify().

In theory, we could go the reverse route: take the descriptor dag and try building a symbolic expression from this which in turn can be fed to Maxima for verification. But I'm far from certain that this is a good idea for those cases which didn't originate in a symbolic expression in the first place. Which brings us back to the request as it currently stands, right?

comment:10 follow-up: Changed 5 years ago by vdelecroix

I just notice that your ticket is based on a very old version of Sage (6.3.beta5)... it would be nice to upgrade to the current 6.5.rc0.

More about the ticket, I am worried because self._symbolic_expression is completely ignored in the following cases

sage: a = QQbar(sqrt(2)) + 1
sage: b = QQbar(sqrt(2)) + QQbar(sqrt(3))

comment:11 in reply to: ↑ 10 ; follow-up: Changed 5 years ago by gagern

Replying to vdelecroix:

More about the ticket, I am worried because self._symbolic_expression is completely ignored in the following cases […]

So what do you propose? That we always try to build a symbolic expression from the descriptor dag?

comment:12 in reply to: ↑ 11 Changed 5 years ago by vdelecroix

Replying to gagern:

Replying to vdelecroix:

More about the ticket, I am worried because self._symbolic_expression is completely ignored in the following cases […]

So what do you propose? That we always try to build a symbolic expression from the descriptor dag?

No. But at least when performing arithmetic operations on QQbar or AA elements the operation are also done on the corresponding symbolic expressions.

comment:13 Changed 5 years ago by git

  • Commit changed from 8597313c855ce7a413cfcabead498059d2cfcbf6 to 998534a33a502a86518f6c541f99cac1ffacdaa5

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

5c81820Drive-by fix to minpoly documentation.
f1bd970Faster exactification using numeric minpoly.
998534aBuild symbolic expression from descriptor graph of algebraic number.

comment:14 follow-up: Changed 5 years ago by gagern

  • Status changed from needs_review to needs_work

I find that it's hard to draw the line. Your example of QQbar(sqrt(2)) + 1 shows that we want symbolic expressions for minpoly construction even if some arguments were not explicitely constructed from symbolic expressions. But where does it stop? In my demo implementation, I decided that caching the symbolic expression after every elementary operation would duplicate a lot of information. So I decided to recursively construct such symbolic expressions. The result is pretty much a construction of symbolic expressions from the descriptor dag, even though I didn't think of it that way when I started coding. It can generate symbolic expression for a lot of things.

And I'm using it unconditionally, which leads to a lot of doctests failing. Some of them for better:

line 3653, rt3b.as_number_field_element():
Expected: … defining polynomial y^4 - 4*y^2 + 1, -a^2 + 2 …
Got:      … defining polynomial y^2 - 3, a …

line 3755, rt2b._exact_value():
Expected: a^3 - 3*a where a^4 - 4*a^2 + 1 = 0 and a in 1.931851652578137?
Got:      a where a^2 - 2 = 0 and a in 1.414213562373095?

line 7171, y.is_complex():
Expected: True
Got:      False

But some of them for worse:

line 360, sage_input(n, verify=True):
Expected: AA(2)
Got:      v = sqrt(AA(2))
          v*v

line 7881, sage_input(n):
Expected: QQbar(sqrt(AA(2))) + 1
Got:      v = sqrt(QQbar(3))
          QQbar(sqrt(AA(2))) + v/v

I don't have any performance comparisons yet, but I guess that there might well be cases where the current approach was fast enough, but the whole detour via symbolic expressions, PARI's algdep and Maxima for verification might have a severe impact.

So my key concern right now is, when do we want to use this approach, and when not to? If multithreading in Python were a more natural concept, I'd be so bold as to suggest that we start both approaches, and the one which terminates first wins, aborting the other. This might make results somewhat unpredictable, though. And given the state of threading support in CPython, I doubt this is realistic in any case.

Here are a few more realistic ways how we could make that case distinctions in a way that might make sense. I'm not completely happy with any of them.

  1. Only do the algdep approach if at least one number in the descriptor dag originated from a symbolic expression
  2. Limit the reconstruction of symbolic expressions to fewer operations, mainly elementary arithmetic operations
  3. Give control to the user: make exactification using algdep not part of the normal exactification procedure, but instead the result of an explicit method invocation
  4. Try to detect exactifications which are likely very slow using the standard approach, e.g. indicated by the expected size of the minpoly for some union field, and try algdep only on those

comment:15 in reply to: ↑ 14 ; follow-up: Changed 5 years ago by vdelecroix

Replying to gagern:

But some of them for worse:

line 360, sage_input(n, verify=True):
Expected: AA(2)
Got:      v = sqrt(AA(2))
          v*v

line 7881, sage_input(n):
Expected: QQbar(sqrt(AA(2))) + 1
Got:      v = sqrt(QQbar(3))
          QQbar(sqrt(AA(2))) + v/v

Not really for worse. The previous doctests assumed that at some point in the code .exactify() would have been called for all nodes of the description tree... Right now, the symbolic code ignore the exactification of the nodes.

EDIT: I removed the non-example I previously wrote here.

Here are a few more realistic ways how we could make that case distinctions in a way that might make sense. I'm not completely happy with any of them.

  1. Only do the algdep approach if at least one number in the descriptor dag originated from a symbolic expression
  2. Limit the reconstruction of symbolic expressions to fewer operations, mainly elementary arithmetic operations
  3. Give control to the user: make exactification using algdep not part of the normal exactification procedure, but instead the result of an explicit method invocation
  4. Try to detect exactifications which are likely very slow using the standard approach, e.g. indicated by the expected size of the minpoly for some union field, and try algdep only on those

Option 3 is definitely needed.

I am against option 1 as we do not want to keep the origin of the objects. By the way, do you have an idea to have quick symbolic representation for

sage: QQbar(sqrt2).sqrt() + QQbar(3).sqrt()

I.e. roots of x^n - a have a simple symbolic counterpart.

For options 2 and 4 there is a right balance to find.

Last edited 5 years ago by vdelecroix (previous) (diff)

comment:16 follow-up: Changed 4 years ago by vdelecroix

Not exactly related to this ticket: I just found out some maxima magic to get the polynomial of symbolic expressions involving square root and such (p. 1041 of the maxima manual):

(%i1) load(to_poly_solve)$
(%i2) first(elim_allbut(first(to_poly(eq,[1,x])),[x]));
                                 4       2
(%o2)                          [x  - 24 x  + 4]
(%i3) eq: x = sqrt(sqrt(3) + sqrt(5) + 1) + sqrt(7);
(%o3)             x = sqrt(7) + sqrt(sqrt(5) + sqrt(3) + 1)
(%i3) first(elim_allbut(first(to_poly(eq, [1,x])), [x]));
         16       14         12          10           8            6
(%o3) [x   - 64 x   + 1648 x   - 23552 x   + 213480 x  - 1191424 x
                                                       4            2
                                            + 3340480 x  - 3562496 x  + 524176]

It is way much much faster than what we have in Sage for minpoly:

sage: maxima.load("to_poly_solve")
"/opt/sage_flatsurf/local/share/maxima/5.35.1/share/to_poly_solve/to_poly_solve.mac"
sage: maxima("eq: x = sqrt(sqrt(3) + sqrt(5) + 1) + sqrt(7);")
x=sqrt(7)+sqrt(sqrt(5)+sqrt(3)+1)
sage: maxima("first(elim_allbut(first(to_poly(eq, [1,x])), [x]));")
[x^16-64*x^14+1648*x^12-23552*x^10+213480*x^8-1191424*x^6+3340480*x^4-3562496*x^2+524176]
sage: %timeit maxima("first(elim_allbut(first(to_poly(eq, [1,x])), [x]));")
10 loops, best of 3: 18 ms per loop
sage: a = sqrt(7) + sqrt(sqrt(5) + sqrt(3) + 1)
sage: %timeit a.minpoly()
1 loops, best of 3: 449 ms per loop

comment:17 in reply to: ↑ 16 Changed 4 years ago by gagern

Replying to vdelecroix:

Not exactly related to this ticket: I just found out some maxima magic to get the polynomial of symbolic expressions involving square root and such

Nice! Do you want to open a spin-off ticket, asking for a way to exploit this functionality in Sage? As you said, it's not exactly on topic here, but it would be a shame if that idea got lost once we have this issue here wrapped up.

comment:18 in reply to: ↑ 15 Changed 4 years ago by gagern

Replying to vdelecroix:

Option 3 is definitely needed.

For options 2 and 4 there is a right balance to find.

I guess once we have a method which can be called by hand, we can use that to experiment a bit, in our daily work, and can ask others to do the same. That might give us some intuition about when to use which. So I'd keep this out of this ticket for now, and instead post a follow-up ticket to investigate that.

comment:19 Changed 4 years ago by git

  • Commit changed from 998534a33a502a86518f6c541f99cac1ffacdaa5 to 96942f09ea8a9cddf2b2dc8fece352935ddfe3b2

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

4b1d01eDrive-by fix to minpoly documentation.
f065ca4Faster exactification using numeric minpoly.
8cc7bc7Build symbolic expression from descriptor graph of algebraic number.
96942f0Trac #16222: Add argument “symbolic” to method “exactify”

comment:20 Changed 4 years ago by gagern

  • Status changed from needs_work to needs_review

Rebased to 6.6.rc2 and added public interface z.exactify(symbolic=True).

I considered “numeric” instead of “symbolic”, but I at least tend to associate “numeric” with “inexact”, so I felt that term might be more confusing. Also considered “algoiritm='symbolic'” but I think in the long run we might have various alternatives which could be tried one after the other, so a flag to enable or disable each sounds more useful. The same holds for a new method as opposed to an argument to exactify. So that's why I chose the interface the way I did.

Now that my new code is not called by default, it shouldn't have any adverse effects on other parts of the code, which might make this easier to review. I hope this can be merged, to enable some evaluation in preparation for #18122 which I filed for automatically choosing between algorithms. I fear that might be hard to get right.

comment:21 Changed 4 years ago by vdelecroix

  • Branch changed from u/gagern/ticket/16222 to u/vdelecroix/16222
  • Commit changed from 96942f09ea8a9cddf2b2dc8fece352935ddfe3b2 to 8773468739544ed9f12b836c7f28a06be37ffd94

Hello,

It's a pity that we got None in

sage: AA(2).sqrt()._descr.quick_symbolic()

I added a new commit to handle that.

On the other hand, did you know that the symbolic ring is very broken to deal with n-th root

sage: a = (-1)^(1/3)
sage: a.n()
0.500000000000000 + 0.866025403784439*I
sage: a**2
1

As bool(I^(2/3) == (-1)^(1/3)) is True we might get into troubles with subtle examples. I am very worried by that.

Vincent


New commits:

8773468Trac 16222: support quick symbolic for roots of x^n -a

comment:22 Changed 4 years ago by vdelecroix

  • Reviewers set to Vincent Delecroix
  • Status changed from needs_review to needs_info

comment:23 follow-up: Changed 4 years ago by vdelecroix

Are there any link/dependency with #14239 or/and #17516?

comment:24 in reply to: ↑ 23 ; follow-up: Changed 4 years ago by gagern

Replying to vdelecroix:

Are there any link/dependency with #14239 or/and #17516?

No, since in general they are not what I'd consider fast. To obtain an exact symbolic number using #14239, we need to have a minimal polynomial, while this one here is aiming to obtain that polynomial more quickly. #17516 would extend #14239 to give results in a larger number of cases, but won't affect this here either. If at all, there is some overlap with #17886, since if that works out relly well, it might make this one here superfluous. I don't think so, though, and since we are only offering this here as an alternative for now, we can use #18122 to evaluate those alternatives once we get there.

… did you know that the symbolic ring is very broken …

No, I did not, but I share your deep concerns. Have you filed a ticket for this yet?

comment:25 in reply to: ↑ 24 Changed 4 years ago by vdelecroix

Replying to gagern:

… did you know that the symbolic ring is very broken …

No, I did not, but I share your deep concerns. Have you filed a ticket for this yet?

#18129

comment:26 Changed 4 years ago by gagern

  • Dependencies set to #15605

I fear we'll have to wait for #15605, the ticket #18129 got duplicated to, before we can trust SR to get this right…

Note: See TracTickets for help on using tickets.