Faster exactification using numeric minpoly

[gagern]

This is a spinoff from #14239 comment:13. 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.

Depends on #15605

Component: number fields

Author: Martin von Gagern

Branch/Commit: u/vdelecroix/16222 @ 8773468

Reviewer: Vincent Delecroix

Issue created by migration from https://trac.sagemath.org/ticket/16222

[gagern]

Branch: u/gagern/ticket/16222

[gagern]

comment:2

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.

2. 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.

3. 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.

4. 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.

5. 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?

6. 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:

bb01e0b	Drive-by fix to minpoly documentation.
491dcc4	Faster exactification using numeric minpoly.

[gagern]

Commit: 491dcc4

[sagetrac-git]

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

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

[sagetrac-git]

Changed commit from 491dcc4 to 8597313

[gagern]

comment:4

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.

[videlec]

comment:6

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

[gagern]

comment:7

Replying to @videlec:

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?

[videlec]

comment:8

Replying to @gagern:

Replying to @videlec:

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().

[gagern]

comment:9

Replying to @videlec:

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?

[videlec]

comment:10

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))

[gagern]

comment:11

Replying to @videlec:

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?

[videlec]

comment:12

Replying to @gagern:

Replying to @videlec:

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.

[sagetrac-git]

Changed commit from 8597313 to 998534a

[sagetrac-git]

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

5c81820	Drive-by fix to minpoly documentation.
f1bd970	Faster exactification using numeric minpoly.
998534a	Build symbolic expression from descriptor graph of algebraic number.

[gagern]

comment:14

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

[videlec]

comment:15

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.

[videlec]

comment:16

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

[gagern]

comment:17

Replying to @videlec:

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.

[gagern]

comment:18

Replying to @videlec:

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.

[sagetrac-git]

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

4b1d01e	Drive-by fix to minpoly documentation.
f065ca4	Faster exactification using numeric minpoly.
8cc7bc7	Build symbolic expression from descriptor graph of algebraic number.
96942f0	Trac #16222: Add argument “symbolic” to method “exactify”

[sagetrac-git]

Changed commit from 998534a to 96942f0

[gagern]

comment:20

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.

[videlec]

Changed commit from 96942f0 to 8773468

[videlec]

comment:21

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:

8773468

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

[videlec]

Changed branch from u/gagern/ticket/16222 to u/vdelecroix/16222

[videlec]

Reviewer: Vincent Delecroix

[videlec]

comment:23

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

[gagern]

comment:24

Replying to @videlec:

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?

[videlec]

comment:25

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

[gagern]

Dependencies: #15605

[gagern]

comment:26

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