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
On Friday, June 11, 2021 at 9:51:46 AM UTC+2 Emmanuel Charpentier wrote:
> 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 <https://cocalc.com> showed the same behavior in
> 9.1 ; so if it is a problem, it is not a recent one…
>
> Is this the expected behavior ?
>
>
--
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sage-support/40475160-972e-4f89-82f7-3e438eed219fn%40googlegroups.com.
