On Fri, 11 Jun 2021, 18:25 Emmanuel Charpentier, <
[email protected]> 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.
---------------------------------------------------------------------------
> 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
> <https://groups.google.com/d/msgid/sage-support/5b389532-91fb-4e33-b8f4-b7ad8e21e419n%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>
--
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/CAAWYfq2CAk82apDHC4AwAacWhvyus1UvFQ2CXVP3px3K5Jgxuw%40mail.gmail.com.
