QQbar operations extremely slow ; is this normal ?

31 views
Skip to first unread message

Emmanuel Charpentier

unread,
Jun 11, 2021, 3:51:46 AMJun 11
to sage-support

As the first part of a demonstration on eigensystems, I was surprised to see that computing QQbar polynomial roots was way slower that computing the roots of the same polynomial expressed as a symbolic expression, or solveing it :

# Relative timings of QQbar.roots and solve
from time import time as stime
BR = QQbar
Dim = 3
set_random_seed(0)
pM = MatrixSpace(BR, Dim).random_element()
sM = matrix([[SR(u.radical_expression()) for u in v] for v in pM])
sl = var("sl")
slI = sl*diagonal_matrix([1]*Dim)
sCM = sM - slI
sCP = sCM.det()
t0 = stime()
SS = solve(sCP, sl)
t1 = stime()
SR = sCP.roots()
t2 = stime()
R1.<pl> = BR[]
plI = pl*diagonal_matrix([1]*Dim)
pCM = pM - plI
pCP = pCM.det()
t3 = stime()
SP = pCP.roots()
t4 = stime()
print("solve : %7.2f, roots(symbolic) : %7.2f, roots(QQbar) : %7.2f"%\
      (t1 -t0, t2 - t1, t4 - t3))

gives :

solve :    2.36, roots(symbolic) :    2.04, roots(QQbar) :  600.07

a quick check on Cocalc showed the same behavior in 9.1 ; so if it is a problem, it is not a recent one…

Is this the expected behavior ?

Emmanuel Charpentier

unread,
Jun 11, 2021, 3:56:14 AMJun 11
to sage-support

Forgot to add : pM.eigenvectors_right() doesn’t returns after 5 hours… And, yes, all pM‘s elements have a radical expression :

sage: sM
[         -sqrt(2) - 1 1/4*I*sqrt(759) - 1/4            -2*sqrt(3)]
[  1/2*I*sqrt(3) + 1/2    1/8*sqrt(33) + 1/8   -1/5*sqrt(29) + 3/5]
[                    0                   1/4  1/2*I*sqrt(23) + 1/2]

HTH,

Fredrik Johansson

unread,
Jun 11, 2021, 9:18:39 AMJun 11
to sage-support
What is reasonable depends on what you expect to be able to do with the solutions. Numerical evaluation should be easy, but if you want canonical forms (minimal polynomials) or the ability to check equalities, that's going to be far more costly.

It looks like the eigenvalues of this matrix will have degree 192 and height several hundred digits. Computing their minimal polynomials explicitly should be doable, but could take minutes, hours or years depending on the algorithm.

The SR solutions just look like the result of expanding the cubic formula; there's not much real computation involved.

The QQbar solutions are huge:

sage: len(str(sage_input(pM.charpoly().roots()[0][0])))
63295

Even so, they are still partially lazy (involving nested field extensions), and reducing them to absolute representations (e.g. with .exactify()) where you can meaningfully test equality is presumably going to take a long time.

How does Calcium fare here? Doing this with arithmetic on minimal polynomials (qqbar_t) is going to take forever. With Calcium fields, representing the solutions lazily is easy. There's no code for representing roots of polynomials implicitly yet, but we can use the cubic formula for a 3x3 matrix. Quick computation with CalciumField in Nemo:

julia> function cubic(C,a,b,c,d)
         a = C(a)
         b = C(b)
         c = C(c)
         d = C(d)
         d0 = b^2 - 3*a*c
         d1 = b*(2*b^2-9*a*c) + 27*a^2*d
         E = ((d1 + sqrt(d1^2 - 4*d0^3))//2) ^ (C(1) // 3)
         D = d0 // E
         w = C(root_of_unity(QQBar, 3))
         w2 = w^2
         x1 = (b + E + D) // (-3*a)
         x2 = (b + w*E + D//w) // (-3*a)
         x3 = (b + w2*E + D//w2) // (-3*a)
         return x1, x2, x3
       end;

julia> function eigproblem()
          C = CalciumField()
          I = C(1im)
          M = C[-sqrt(C(2)) - 1  C(1)//4*I*sqrt(C(759)) - C(1)//4  -2*sqrt(C(3));
                 C(1)//2*I*sqrt(C(3)) + C(1)//2  C(1)//8*sqrt(C(33)) + C(1)//8  -C(1)//5*sqrt(C(29)) + C(3)//5;
                 C(0)  C(1)//4  C(1)//2*I*sqrt(C(23)) + C(1)//2]
          CP = charpoly(PolynomialRing(C, "x")[1], M)
          return CP, cubic(C, coeff(CP, 3), coeff(CP, 2), coeff(CP, 1), coeff(CP, 0));
       end;

julia> @time CP, (x1, x2, x3) = eigproblem();
  0.215785 seconds (135.20 k allocations: 11.313 MiB, 46.95% gc time)

julia> x1
0.108407 + 1.67305*I {(-160*a^2+20*a*d+80*a*f*i-160*a*h-60*a+60*d*f*g-50*d*f*i-20*d*h-15*d+24*e-80*f*h*i-150*f*i+60*g*i-420*h+213)/(480*a) where a = 3.36525 - 3.03539*I [Pow(-54.9070 - 75.1596*I {(6*b+90*d*e-1800*d*f*g*h-3375*d*f*g+3000*d*f*h*i+4950*d*f*i-41175*d*g*i+2250*d*h-42045*d+360*e*f*i+1440*e*h+1890*e+9225*f*g+3600*f*h*i+5270*f*i-1800*g*h*i-3375*g*i+10800*g+40930*h+32400*i+71955)/3200}, 0.333333 {1/3})], b = 13494.0 - 27535.7*I [Sqrt(-5.76129e+8 - 7.43136e+8*I {-711000*d*e*f*g*h+385050*d*e*f*g+1026000*d*e*f*h*i+241200*d*e*f*i-5247000*d*e*g*h*i-7971750*d*e*g*i+54000*d*e*g-5202300*d*e*h+162000*d*e*i-7560900*d*e-1440000*d*f*g*h*i-13943250*d*f*g*h-3105000*d*f*g*i-23751150*d*f*g+22459500*d*f*h*i-8640000*d*f*h-113185725*d*f*i-14985000*d*f-38975250*d*g*h*i+1350000*d*g*h-107574750*d*g*i+48888000*d*g+4050000*d*h*i-197409600*d*h-149796000*d*i-341006175*d+129000*e*f*g*h+216000*e*f*g*i+272250*e*f*g+2902800*e*f*h*i+3985200*e*f*i-648000*e*f-711000*e*g*h*i+864000*e*g*h-25834950*e*g*i+1134000*e*g+2592000*e*h*i+5213700*e*h+3402000*e*i+34315260*e+2160000*f*g*h*i+297276750*f*g*h+19767000*f*g*i+524963250*f*g+119813100*f*h*i-6480000*f*h+238612275*f*i+7119000*f-475643250*g*h*i+27798000*g*h+2049398850*g*i+49248000*g+70434000*h*i-343497600*h+123444000*i-1837827855})], c = 27.5500 [c^2-759=0], d = 5.74456 [d^2-33=0], e = 5.38516 [e^2-29=0], f = 4.79583 [f^2-23=0], g = 1.73205 [g^2-3=0], h = 1.41421 [h^2-2=0], i = I [i^2+1=0]}

julia> CP(x1)
0e-21 - 0e-21*I {(-1024000*a^6+57600*a^3*d*e-1152000*a^3*d*f*g*h-2160000*a^3*d*f*g+1920000*a^3*d*f*h*i+3168000*a^3*d*f*i-26352000*a^3*d*g*i+1440000*a^3*d*h-26908800*a^3*d+230400*a^3*e*f*i+921600*a^3*e*h+1209600*a^3*e+5904000*a^3*f*g+2304000*a^3*f*h*i+3372800*a^3*f*i-1152000*a^3*g*h*i-2160000*a^3*g*i+6912000*a^3*g+26195200*a^3*h+20736000*a^3*i+46051200*a^3-907200*d*e*f*g*h+568080*d*e*f*g+907200*d*e*f*h*i+201600*d*e*f*i-4017600*d*e*g*h*i-7484400*d*e*g*i-3238560*d*e*h-5720220*d*e-34268400*d*f*g*h+21208410*d*f*g+17771400*d*f*h*i-269751300*d*f*i+155350800*d*g*h*i+109899450*d*g*i+28690380*d*h+72015435*d-1252800*e*f*g*h-745200*e*f*g+2842560*e*f*h*i+2160000*e*f*i-907200*e*g*h*i-81511920*e*g*i-12800160*e*h+106462116*e-274287600*f*g*h-477318150*f*g-365157880*f*h*i-728050500*f*i+1689411600*g*h*i-88231590*g*i-1638440820*h+2202115707)/(27648000*a^3) where a = 3.36525 - 3.03539*I [Pow(-54.9070 - 75.1596*I {(6*b+90*d*e-1800*d*f*g*h-3375*d*f*g+3000*d*f*h*i+4950*d*f*i-41175*d*g*i+2250*d*h-42045*d+360*e*f*i+1440*e*h+1890*e+9225*f*g+3600*f*h*i+5270*f*i-1800*g*h*i-3375*g*i+10800*g+40930*h+32400*i+71955)/3200}, 0.333333 {1/3})], b = 13494.0 - 27535.7*I [Sqrt(-5.76129e+8 - 7.43136e+8*I {-711000*d*e*f*g*h+385050*d*e*f*g+1026000*d*e*f*h*i+241200*d*e*f*i-5247000*d*e*g*h*i-7971750*d*e*g*i+54000*d*e*g-5202300*d*e*h+162000*d*e*i-7560900*d*e-1440000*d*f*g*h*i-13943250*d*f*g*h-3105000*d*f*g*i-23751150*d*f*g+22459500*d*f*h*i-8640000*d*f*h-113185725*d*f*i-14985000*d*f-38975250*d*g*h*i+1350000*d*g*h-107574750*d*g*i+48888000*d*g+4050000*d*h*i-197409600*d*h-149796000*d*i-341006175*d+129000*e*f*g*h+216000*e*f*g*i+272250*e*f*g+2902800*e*f*h*i+3985200*e*f*i-648000*e*f-711000*e*g*h*i+864000*e*g*h-25834950*e*g*i+1134000*e*g+2592000*e*h*i+5213700*e*h+3402000*e*i+34315260*e+2160000*f*g*h*i+297276750*f*g*h+19767000*f*g*i+524963250*f*g+119813100*f*h*i-6480000*f*h+238612275*f*i+7119000*f-475643250*g*h*i+27798000*g*h+2049398850*g*i+49248000*g+70434000*h*i-343497600*h+123444000*i-1837827855})], c = 27.5500 [c^2-759=0], d = 5.74456 [d^2-33=0], e = 5.38516 [e^2-29=0], f = 4.79583 [f^2-23=0], g = 1.73205 [g^2-3=0], h = 1.41421 [h^2-2=0], i = I [i^2+1=0]}

These representations are pretty similar to the SR solutions, just displayed more compactly. You can use them for numerics:

julia> AcbField(1000)(CP(x1))
[+/- 4.69e-648] + [+/- 4.25e-648]*im

julia> AcbField(10000)(CP(x1))
[+/- 7.21e-6336] + [+/- 7.18e-6336]*im


The following exact test fails because the absolute representations of the algebraic numbers are too large:

julia> CP(x1) == 0
ERROR: Unable to perform operation (failed deciding truth of a predicate): isequal


It's an interesting problem to make this work (in reasonable time).

Fredrik

Emmanuel Charpentier

unread,
Jun 11, 2021, 1:25:43 PMJun 11
to sage-support
OK. That was not an error of mine (nor my Sage installation), but a a genuine problem. So no indication to file a ticket.

Note : Numerics with"standard" are not a problem :

sage: %time matrix([[CDF(u) for u in v] for v in pM]).eigenvectors_right()
CPU times: user 1.95 ms, sys: 15 µs, total: 1.96 ms
Wall time: 1.85 ms
[(-1.5855741097050124 - 1.9835895314317984*I,
  [(0.9524479486112228, -0.2837102376875301 - 0.11002559239003057*I, 0.011400252227268036 - 0.010761481586469179*I)],
  1),
 (0.10840730449357763 + 1.6730496133213792*I,
  [(0.8835538411020463, 0.2631022513327894 - 0.3627905504445216*I, 0.05890961187635257 + 0.12256626515553348*I)],
  1),
 (0.4060235736555931 + 2.708455679766778*I,
  [(0.8864942697472706, 0.26813150779882816 - 0.24992379361655742*I, -0.24416421274380526 - 0.14196949964794017*I)],
  1)]

But if you need extended precision, there's a snag :

```
sage: C= ComplexField(200)
sage: %time foo = matrix([[C(u) for u in v] for v in pM]).eigenvectors_right()
<timed exec>:1: UserWarning: Using generic algorithm for an inexact ring, which may result in garbage from numerical precision issues.
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<timed exec> in <module>

/usr/local/sage-9/local/lib/python3.9/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.eigenvectors_right (build/cythonized/sage/matrix/matrix2.c:47695)()
   6709             implemented for RDF and CDF, but not for Rational Field
   6710         """
-> 6711         return self.transpose().eigenvectors_left(other=other, extend=extend)
   6712 
   6713     right_eigenvectors = eigenvectors_right

/usr/local/sage-9/local/lib/python3.9/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.eigenvectors_left (build/cythonized/sage/matrix/matrix2.c:46717)()
   6623         from sage.rings.qqbar import QQbar
   6624         from sage.categories.homset import hom
-> 6625         eigenspaces = self.eigenspaces_left(format='galois', algebraic_multiplicity=True)
   6626         evec_list=[]
   6627         n = self._nrows

/usr/local/sage-9/local/lib/python3.9/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.eigenspaces_left (build/cythonized/sage/matrix/matrix2.c:42000)()
   6127             msg = ("eigenspaces cannot be computed reliably for inexact rings such as {0},\n",
   6128                    "consult numerical or symbolic matrix classes for other options")
-> 6129             raise NotImplementedError(''.join(msg).format(self.base_ring()))
   6130 
   6131         format = self._eigenspace_format(format)

NotImplementedError: eigenspaces cannot be computed reliably for inexact rings such as Complex Field with 200 bits of precision,
consult numerical or symbolic matrix classes for other options
```

so you're back to "manual" computation of eigenvalues, then eigenvectors, and maintaining precision sin't obvious...

Thank you very much.

Dima Pasechnik

unread,
Jun 11, 2021, 2:22:43 PMJun 11
to sage-support


On Fri, 11 Jun 2021, 18:25 Emmanuel Charpentier, <emanuel.c...@gmail.com> wrote:
OK. That was not an error of mine (nor my Sage installation), but a a genuine problem. So no indication to file a ticket.

Note : Numerics with"standard" are not a problem :

sage: %time matrix([[CDF(u) for u in v] for v in pM]).eigenvectors_right()
CPU times: user 1.95 ms, sys: 15 µs, total: 1.96 ms
Wall time: 1.85 ms
[(-1.5855741097050124 - 1.9835895314317984*I,
  [(0.9524479486112228, -0.2837102376875301 - 0.11002559239003057*I, 0.011400252227268036 - 0.010761481586469179*I)],
  1),
 (0.10840730449357763 + 1.6730496133213792*I,
  [(0.8835538411020463, 0.2631022513327894 - 0.3627905504445216*I, 0.05890961187635257 + 0.12256626515553348*I)],
  1),
 (0.4060235736555931 + 2.708455679766778*I,
  [(0.8864942697472706, 0.26813150779882816 - 0.24992379361655742*I, -0.24416421274380526 - 0.14196949964794017*I)],
  1)]

But if you need extended precision, there's a snag :

```
sage: C= ComplexField(200)
sage: %time foo = matrix([[C(u) for u in v] for v in pM]).eigenvectors_right()
<timed exec>:1: UserWarning: Using generic algorithm for an inexact ring, which may result in garbage from numerical precision issues.

in fact, arb has some linear algebra implemented, only that functionality is not wrapped in Sage.


--
You received this message because you are subscribed to the Google Groups "sage-support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sage-support...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sage-support/5b389532-91fb-4e33-b8f4-b7ad8e21e419n%40googlegroups.com.

Vincent Delecroix

unread,
Jun 11, 2021, 3:27:05 PMJun 11
to sage-s...@googlegroups.com
hmmm, 10 years ago it was indeed not in sage but

sage: sM = matrix(3, [-sqrt(2) - 1, 1/4*I*sqrt(759) - 1/4, -2*sqrt(3),
....: 1/2*I*sqrt(3) + 1/2, 1/8*sqrt(33) + 1/8, -1/5*sqrt(29) + 3/5,
....: 0, 1/4, 1/2*I*sqrt(23) + 1/2])
sage: sM.change_ring(ComplexBallField(256)).eigenvectors_right_approx()
<ipython-input-2-597c50ed1766>:1: FutureWarning: This
class/method/function is marked as experimental. It, its functionality
or its interface might change without a formal deprecation.
See https://trac.sagemath.org/30393 for details.
sM.change_ring(ComplexBallField(Integer(256))).eigenvectors_right_approx()
[([-1.585574109705012818320754555380077143033164486087030751802791286603000260637687
+/- 3.06e-79] -
[1.983589531431798407272840921819657822844388213332549281429493414447825431290207
+/- 3.05e-79]*I,

[([-0.7645889954375910796175588680760753644497544472611288268496784517551156034840851
+/- 3.49e-80] +
[0.5679443307838031693126382504352883182783867453515207101490113908352846374511589
+/- 4.71e-80]*I,
[0.2933600072059874975695356608056928703375239897336479003693686118090527352449637
+/- 2.62e-80] -
[0.08085193950422595460559726747128398322364244288273211664977438986727139956036997
+/- 1.86e-81]*I,
[-0.002734621817507568627793755191059520745414855464984622057215624183928681966128039
+/- 3.48e-83] +
[0.01543687404549433149004830726668865175938723340836903249472149019570000983681407
+/- 6.01e-82]*I)],
1),

([0.1084073044935780329693711992135459902927450863702466775904749191406132984593997
+/- 2.35e-80] +
[1.673049613321379085900796686139959744878064234423412272818114017801709044221350
+/- 4.95e-79]*I,

[([-0.2673371902772041016297440268511643741966917470300012123253823050737523650735768
+/- 3.44e-80] -
[0.02255276973433239605480951065284270406141140134762632616605864073042511737197782
+/- 6.64e-82]*I,
[-0.08886719147182478198522558738553193426651781904829084677832088009136722752956601
+/- 2.61e-81] +
[0.1030539596890842335569665014655476103537849006257746372204680855562700952313454
+/- 1.45e-80]*I,
[-0.01469578961696251938128853753478067963226607548513094281484268115068461631364436
+/- 3.99e-82] -
[0.03858858880485532751435734250502627274540851776921042689105965744221705829132246
+/- 3.42e-81]*I)],
1),

([0.4060235736555933190310210654841992389482805815876850469981213132923920738981516
+/- 2.51e-80] +
[2.708455679766779092170763267761045037964677499861201681854033953870245005022438
+/- 1.50e-79]*I,

[([0.2103718533466063742076462260260793168410401002165862474556218965077763660977134
+/- 3.27e-80] +
[0.05448236953479104093710977337474251123511802942923003190706717417804356762302697
+/- 3.17e-81]*I,
[0.07898952661654168047766696024137805415491440866540645841786257582178406891492673
+/- 1.47e-81] -
[0.04282993479195057509885682538044990967004220531263985306993932362186424999425426
+/- 2.1e-84]*I,
[-0.04921683614015854406786925932481355184588924524472501364264944877387240329246417
+/- 4.04e-82] -
[0.04869634593105051448079497440445424999631243725298577782757982000687160120825688
+/- 3.53e-81]*I)],
1)]

Emmanuel Charpentier

unread,
Jun 11, 2021, 4:55:21 PMJun 11
to sage-support

Very nice ! I wasn’t aware of it… Now a couple questions/remarks :

  • Is there a reason for this being implemented for ComplexBallField but not for ComplexField (except CDF) ? What about ComplexIntervalField ?

  • I am not aware of anything comparing the abilities of the various numerical tools available in Sagemath in the documentation. Did I miss it ?

At the very least, a tutorial illustrating them would be useful. A better (but possibly quite heavier) text would be additions to the numerical computing chapters of a new edition of Computational mathematics with Sagemath, which is more and more in order…

What do you think ?

Vincent Delecroix

unread,
Jun 11, 2021, 5:27:19 PMJun 11
to sage-s...@googlegroups.com
You can ask Thierry to finish his ticket

https://trac.sagemath.org/ticket/15944

It was started 7 years ago and the description says "The doc is not
yet finished, but should be quite soon."

That would be a good start :)

Le 11/06/2021 à 22:55, Emmanuel Charpentier a écrit :
>
>
> Very nice ! I wasn’t aware of it… Now a couple questions/remarks :
>
> -
>
> Is there a reason for this being implemented for ComplexBallField but
> not for ComplexField (except CDF) ? What about ComplexIntervalField ?
> -
>
> I am not aware of anything comparing the abilities of the various
> numerical tools available in Sagemath in the documentation. Did I miss it ?
>
> At the very least, a tutorial illustrating them would be useful. A better
> (but possibly quite heavier) text would be additions to the numerical
> computing chapters of a new edition of *Computational mathematics with
> Sagemath*, which is more and more in order…
Reply all
Reply to author
Forward
0 new messages