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.
Le vendredi 11 juin 2021 à 15:18:39 UTC+2, Fredrik Johansson a écrit :
> 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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/sage-support/5b389532-91fb-4e33-b8f4-b7ad8e21e419n%40googlegroups.com.
