QQbar benchmarks

[Fredrik Johansson]

Hi all,

I'm looking for benchmark problems for QQbar/AA arithmetic. Ideally such a problem will:

* Be reducible to a short program that, apart from using QQbar/AA operations, is reasonably self-contained.
* Reflect real-world use, i.e. originate from using Sage to solve an actual mathematical problem, with QQbar/AA arithmetic being a significant ingredient.

I'm interested in such benchmark problems to evaluate how well Calcium (https://fredrikj.net/calcium/) would work to implement the field of algebraic numbers. As such, I'd like to identify any situations where Calcium is slower than Sage's QQbar or where functionality missing.

There are actually two independent types that may serve to implement the field of algebraic numbers: qqbar_t and ca_t. A Calcium qqbar_t uses an "absolute" representation (minimal polynomial), which is nice since results are kept reduced and canonical, but can be slow for field arithmetic. A Calcium ca_t uses a "relative" representation involving symbolic expressions and number fields, somewhat similar to Sage's QQbar or Magma's algebraically closed field. The ca_t type should give the best performance in most cases, but qqbar_t may be better in some circumstances.

Some benchmark problems I've collected so far:

huge_expr (https://ask.sagemath.org/question/52653/)
====================================================
sage qqbar:         ~forever (hours)
calcium qqbar:      8.5 s
calcium ca:         8.1 s

heptadecagon (sage/rings/qqbar.py)
====================================================
sage qqbar:         ~forever (hours)
calcium qqbar:      2.11 s
calcium ca:         0.011 s

heptadecagon2 (sage/rings/qqbar.py)
====================================================
sage qqbar:         0.060 s
calcium qqbar:      2.11 s
calcium ca:         0.0031 s

symbench_R1 (https://wiki.sagemath.org/symbench)
====================================================
sage qqbar (numerical answer):  0.0041 s
sage qqbar (exactify):          0.049 s

calcium qqbar:                  0.095 s

calcium ca (numerical answer):  0.0010 s

calcium ca (qqbar output)       0.057 s

Fredrik

[Vincent Delecroix]

Dear Fredrik,

One technical question: I thought that your ca_t implementation
used multivariate polynomials. This is what Magma does but not
what sage does. The latter uses expression trees and take union
fields anytime there is an exactification to be done. Is it a
misconception of mine?


For "real world" problems

1) algebraic solutions of zero dimensional varieties
defined over QQ

sage: R.<x,y,z> = QQ[]
sage: A = x^4 + y^4 + z^3 - 1
sage: B = x^2*z + 2*y^2*x - x*y*z + 5
sage: C = y*z^2 - 2*(x+y+z) + 2
sage: %time solutions = R.ideal([A,B,C]).variety(QQbar)
CPU times: user 24.3 s, sys: 23.3 ms, total: 24.3 s
Wall time: 24.4 s

Though here is the bottleneck is Groebner basis computation
of multivariate polynomials and not really QQbar (21s over the 24s).


2) Two examples somehow similar to your huge_expr

sage: a = sum(AA(p).sqrt() for p in prime_range(18))
sage: %time a.exactify()
CPU times: user 30.7 s, sys: 13.3 ms, total: 30.8 s
Wall time: 30.8 s
sage: a.degree() # does not terminates

In the above situation, the problem is that sage is representing
the number a as a complicated element in a rather complicated
number field. Computing the minpoly involves apparently too
complicated linear algebra (?). However, if instead of exactifying
in that way but using the minimal polynomial (which looks reasonable
in this range of degrees) the latter would be very quick.

sage: b = sum(QQbar.zeta(n) for n in range(1,8))
sage: b.exactify() # does not terminate

Here sage is stuck at building union fields.


Would be nice to have calcium used in sage!

Best
Vincent

[Fredrik Johansson]

Hi Vincent,

On Sunday, April 25, 2021 at 4:26:18 PM UTC+2 vdelecroix wrote:

Dear Fredrik,

One technical question: I thought that your ca_t implementation
used multivariate polynomials. This is what Magma does but not
what sage does. The latter uses expression trees and take union
fields anytime there is an exactification to be done. Is it a
misconception of mine?

Yes, ca_t uses multivariate polynomials (in fact, multivariate fraction field elements). It uses univariate polynomials in the special case of simple (univariate) number fields.

Symbolic expressions are effectively used to represent non-field operations generating new extension elements (e.g. square roots, when not converted to absolute qqbar elements).

With qqbar_t:

In [2]: %time sum(qqbar(p).sqrt() for p in [2, 3, 5, 7, 11, 13, 17])
CPU times: user 286 ms, sys: 4.23 ms, total: 291 ms
Wall time: 290 ms
Out[2]: 19.0734 (deg 128)

In [3]: %time sum(qqbar(RootOfUnity(n)) for n in range(1,8))
CPU times: user 24.1 ms, sys: 30 µs, total: 24.1 ms
Wall time: 23.5 ms
Out[3]: 0.932507 + 4.46494*I (deg 96)

With ca_t:

In [2]: %time sum(ca(p).sqrt() for p in [2, 3, 5, 7, 11, 13, 17])
CPU times: user 2.77 ms, sys: 265 µs, total: 3.03 ms
Wall time: 2.63 ms
Out[2]: 19.0734 {a+b+c+d+e+f+g where a = 4.12311 [a^2-17=0], b = 3.60555 [b^2-13=0], c = 3.31662 [c^2-11=0], d = 2.64575 [d^2-7=0], e = 2.23607 [e^2-5=0], f = 1.73205 [f^2-3=0], g = 1.41421 [g^2-2=0]}

In [3]: %time sum(ca(qqbar(RootOfUnity(n))) for n in range(1,8))
CPU times: user 3.13 ms, sys: 0 ns, total: 3.13 ms
Wall time: 2.55 ms
Out[3]: 0.932507 + 4.46494*I {a+b+c+d where a = 0.623490 + 0.781831*I [a^6+a^5+a^4+a^3+a^2+a+1=0], b = 0.309017 + 0.951057*I [b^4+b^3+b^2+b+1=0], c = 1.73205*I [c^2+3=0], d = I [d^2+1=0]}

Of course, it is debatable whether the result has really been "computed" here. If you convert the result to a qqbar_t, it will basically take as long as using qqbar_t in the first place. But these examples are trivial; in other cases, the multivariate arithmetic will be helpful.

Fredrik

[Michael Orlitzky]

On Sun, 2021-04-25 at 15:07 +0200, Fredrik Johansson wrote:
> Hi all,
>
> I'm looking for benchmark problems for QQbar/AA arithmetic. Ideally such a
> problem will:
>
> * Be reducible to a short program that, apart from using QQbar/AA
> operations, is reasonably self-contained.
> * Reflect real-world use, i.e. originate from using Sage to solve an actual
> mathematical problem, with QQbar/AA arithmetic being a significant
> ingredient.

I'm interested in Gram-Schmidt and the coordinate_vector() methods that
solve linear systems over AA.

For a benchmark, you could conjure a set of long vectors with entries
in AA/QQbar, and orthonormalize them using the kindergarten algorithm
that everyone knows. That will leave a linearly-independent set, and if
you keep the input/output indices matched up, will let you know which
of the input vectors correspond to that linearly independent set. Call
the set of those input vectors "X". Then to benchmark the
coordinate_vector() methods, you can generate random elements in
W=span(X), and solve the system to find their W-coordinates in terms of
the basis X. (Using the de-orthonormalized vectors as the basis makes
the system harder to solve, and more representative of the general
case.)


Unimportant background information:

A Euclidean Jordan algebra (EJA) is a real vector space where every
element "x" has a self-adjoint "left multiplication by x" operation
that I'll write L_x.

It turns out that every EJA is a mishmash of various Hermitian matrix
algebras, which is nice, because it lets you think of their elements as
being some concrete object rather than just an abstract vector space
element. For example, the space of real symmetric matrices with

L_x(y) := (xy + yx)/2

induces an EJA. But when people think of a symmetric matrix, they're
thinking of it in terms of the standard basis, which leads to a
problem; namely that the basis

[1 0] [0 1] [0 0]
[0 0], [1 0], [0 1]

for S^2 is not normalized. This raises a problem: while L_x is always
self-adjoint, its matrix representation won't be symmetric unless your
basis is orthonormal. And you need its matrix to be symmetric if you
want to do anything cool with it in sage, because sage assumes that
when you have a matrix full of numbers, they're with respect to the
standard basis.

In short, this means that a sqrt(2) must become involved, and the "2"
isn't special. Here's where Gram-Schmidt comes in. Afterwards, matrices
aren't vectors in sage, so converting back and forth between
matrix/vector representations involves a lot of to_vector() and
coordinate_vector() over AA.

Things get worse: when constructing a subalgebra, its basis needs to be
orthonormalized as well, and we need to be able to convert back and
forth between the orthonormalized superalgebra coordinates and the
orthonormalized subalgebra coordinates. Thus we wind up
orthonormalizing vectors that already have ugly entries, and solving
systems with uglier and uglier entries in AA. This is a bottleneck when
the algebras get "big," which they do even in simple examples (the
Albert algebra has dimension 27).

[John Cremona]

+1 for Vincent's zero-dimensional variety point finding. I used to
use Magma's AlgebraicallyClosedField a lot for this (for some project
which I did before Sage existed). The idea is to find the points over
Qbar, then construct the absolute field they generate, then find them
again over that field before continuing. I can imagine a nice Sage
function to do that automatically including a field_of_definition_of
points method on such a variety.

John

> --
> You received this message because you are subscribed to the Google Groups "sage-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/1620f0b7f268a16e477d74f7295d9b1340285604.camel%40orlitzky.com.

[Fredrik Johansson]

On Sunday, April 25, 2021 at 4:26:18 PM UTC+2 vdelecroix wrote:

Dear Fredrik,

One technical question: I thought that your ca_t implementation
used multivariate polynomials. This is what Magma does but not
what sage does. The latter uses expression trees and take union
fields anytime there is an exactification to be done. Is it a
misconception of mine?


For "real world" problems

1) algebraic solutions of zero dimensional varieties
defined over QQ

sage: R.<x,y,z> = QQ[]
sage: A = x^4 + y^4 + z^3 - 1
sage: B = x^2*z + 2*y^2*x - x*y*z + 5
sage: C = y*z^2 - 2*(x+y+z) + 2
sage: %time solutions = R.ideal([A,B,C]).variety(QQbar)
CPU times: user 24.3 s, sys: 23.3 ms, total: 24.3 s
Wall time: 24.4 s

Though here is the bottleneck is Groebner basis computation
of multivariate polynomials and not really QQbar (21s over the 24s).

I managed to isolate the final QQbar part of this computation in a way that I could benchmark outside of Sage. Results:

sage qqbar: 500 ms
calcium ca: 80 ms

calcium qqbar: ~forever (hours)

What takes forever here is evaluating polynomials with coefficients in Q at QQbar arguments, using a generic algorithm (Horner's rule or adding powers). This is slow in the minpoly representation, because you get huge temporary results which each need to be factored.  You could throw in some specialized polynomial evaluation code that generates a number field element and then computes the absolute minimal polynomial, but it'd essentially be reinventing the number field arithmetic which Sage's QQbar and Calcium's ca do automatically.

I'd be interested in some way to improve arithmetic in the minpoly representation when the operands lie in the same field (e.g. by quickly recognizing such operands using LLL). My earlier attempts to implement such hacks were unsuccessful (only slowed things down on examples I tried), but perhaps it can make difference here if I get it right.

Fredrik

[Samuel Lelievre]

Hi Fredrik,

These papers might provide benchmark problems

for QQbar/AA arithmetic which originate in actual

math problems, though not in Sage.

Pierre-Vincent Koseleff, Fabrice Rouillier, Cuong Tran
On the Sign of a Trigonometric Expression
ISSAC’15, July 6–9, 2015, Bath, United Kingdom.
http://dx.doi.org/10.1145/2755996.2756664.

Koseleff, P.-V., Pecker, D., and Rouillier, F.
Computing Chebyshev knot diagrams.
http://arxiv.org/abs/1001.5192

Koseleff, P.-V., Pecker, D., and Rouillier, F.

The first rational Chebyshev knots.

Journal of Symbolic Computation 45 (2010), 1341–1358.

See Pierre-Vincent Koseleff's publications page:

https://webusers.imj-prg.fr/~pierre-vincent.koseleff/publications.html

There are some interesting algebraic numbers in Theorem 3 of

Shigeki Akiyama

Minimum polyhedron with n vertices

https://arxiv.org/abs/1704.03149

which maybe can also provide benchmark problems.

Happy benchmarking!   --Samuel