[sage-support] Re: Function defined with matrices won't plug in values for function
Hi Saad, sorry that it took long to answer; I thought others would reply sooner than I. On 2018-12-05, saad khalid wrote: > If I run it, it gives me results with the error: " > > UserWarning: Using generic algorithm for an inexact ring, which will probably > give incorrect results due to numerical precision issues." It's not an error but a warning. The absence of a warning doesn't mean that the result is more trustworthy than a result of a computation that doesn't create a warning. See below: When doing numerical computations, it is advisable to assess the results by cross verifications and consistency checks. > I expected this error to be some minute difference in the far out decimal > place range, like I so often see with > complex numbers, but the answer I get is actually Quite a bit off from the > correct answer. No surprise. That's in the nature of numerical computations. > My question is, why is the default an alg that seems to give incorrect > answers, instead of just defaulting to CDF? > Also, is there a better way of doing this besides just appending CDF to the > beginning every time? The default complex field CC has a precision of 53 bits. In principle you can create a complex field with arbitrary finite precision, by using ComplexField(bitwise_precision) CDF is just a wrapper of machine double precision complex numbers, which have a fixed precision of 53 bit as well. However, computations in CDF are based on the GNU Scientific Library (GSL). Apparently CC uses a generic algorithm to compute eigenvalues --- which is both slower and more prone to rounding effects than the algorithm in GSL. So, as usual when doing numerical computations, one needs to make the best of the available resources. For example, if in some parts of your computation you do "ordinary" arithmetic computations and need to do it in very high precision, you could create a complex field in the required precision, and build your matrix M on top of it; and when you in another part of your computations you need a stable numerical algorithm that isn't offered by CC, then you can always change to CDF. For example: sage: C = ComplexField(500) sage: M = random_matrix(C, 10) sage: M2 = M.change_ring(CDF) sage: abs(M.determinant()^20 - (M^20).determinant()) 3.22553502188764194905953565908366462673376034391868764719927365055487394028435964347935057202262954253087699368051026966711621636866859093419028355926e-59 sage: abs(M2.determinant()^20 - (M2^20).determinant()) 2.1231804490158553e+41 Of course, the absolute value of the difference should theoretically be zero. So, as you can see here, in my example it is many orders of magnitude better to use a high precision field instead of CDF. Also, the eigenvalues and characteristic polynomial are more consistent than those computed in CDF: Let's compare the characteristic polynomial with the product of (x-eigenvalue): sage: R. = C[] sage: max(abs(c) for c in (M.charpoly()-prod((x-ev) for ev in M.eigenvalues())).list()) /home/king/Sage/git/sage/src/bin/sage-ipython:1: UserWarning: Using generic algorithm for an inexact ring, which will probably give incorrect results due to numerical precision issues. #!/usr/bin/env sage-python23 4.69238225433539279163184034797611678687652260951909375574059096465606393834197062303005271312591139060256023453735495755056822286510738017751316446341e-148 sage: max(abs(c) for c in (M2.charpoly()-prod((x-ev) for ev in M2.eigenvalues())).list()) 3.096434044909977e-12 But if the precision of the field C is only marginally larger than the precision of CDF, the picture changes: sage: C = ComplexField(60) sage: R. = C[] sage: M = random_matrix(C, 10) sage: M2 = M.change_ring(CDF) sage: abs(M.determinant()^20 - (M^20).determinant()) 7.9015802853808835e72 sage: abs(M2.determinant()^20 - (M2^20).determinant()) 4.509812763227185e+43 sage: max(abs(c) for c in (M.charpoly()-prod((x-ev) for ev in M.eigenvalues())).list()) /home/king/Sage/git/sage/src/bin/sage-ipython:1: UserWarning: Using generic algorithm for an inexact ring, which will probably give incorrect results due to numerical precision issues. #!/usr/bin/env sage-python23 9.9301366129890920e-16 sage: max(abs(c) for c in (M2.charpoly()-prod((x-ev) for ev in M2.eigenvalues())).list()) 3.5482789005657048e-12 This time, determinant and matrix multiplication are more consistent in CDF than in C (although the precision used in C is still slightly better than in CDF!). But the consistency of eigenvalue and charpoly computation is still better in C and in CDF. However, the fact that eigenvalues and charpoly are more *consistent* in CC than in CDF does not mean that the eigenvalues are more close to reality. So, let's start with a matrix with known eigenvalues, then compare with computations in CC resp. in CDF resp. in precision 2000: sage: EV = [CC.random_element() for _ in range(10)] sage:
[sage-support] Re: Function defined with matrices won't plug in values for function
On 2018-12-15, Simon King wrote: > It's not an error but a warning. The absence of a warning doesn't mean > that the result is more trustworthy than a result of a computation that > doesn't create a warning. Ooops, one negation too many. I meant to say: "The absence of a warning doesn't mean that the result is more trustworthy than the result of a computation that DOES create a warning. Sorry. -- 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 post to this group, send email to sage-support@googlegroups.com. Visit this group at https://groups.google.com/group/sage-support. For more options, visit https://groups.google.com/d/optout.
[sage-support] Re: Function defined with matrices won't plug in values for function
On Wednesday, December 5, 2018 at 10:15:52 AM UTC-8, saad khalid wrote:
>
> Thank you all for the thoughtful responses, your comments helped me a lot.
> One last bit of trouble that I've been having related to this topic. This
> was is the code I was running:
>
> Mz = matrix([[0,1],[1,0]])
> Mx = matrix([[1,0],[0,-1]])
> M1 = Matrix([[1,0],[0,1]])
> import numpy as numpy
> #R. = QQ[]
> _ = var('J,hx,alpha')
> alphas = numpy.linspace(0,1,10)
> Bp = Mz.tensor_product(M1).tensor_product(M1).tensor_product(M1) * M1.
> tensor_product(Mz).tensor_product(M1).tensor_product(M1) * M1.
> tensor_product(M1).tensor_product(Mz).tensor_product(M1) * M1.
> tensor_product(M1).tensor_product(M1).tensor_product(Mz)
> Plaqz = Mx.tensor_product(M1).tensor_product(M1).tensor_product(M1)
> F = lambda s: -(1-s)*Bp - s*Plaqz
>
> show(F(.1).eigenvalues())
>
> If I run it, it gives me results with the error: "
>
> UserWarning: Using generic algorithm for an inexact ring, which will probably
> give incorrect results due to numerical precision issues."
>
> I expected this error to be some minute difference in the far out decimal
> place range, like I so often see with
> complex numbers, but the answer I get is actually Quite a bit off from the
> correct answer. Upon googling, i found that this can
> be fixed by appending CDF to the beginning of every matrix definition, like
> this:
>
>
> My question is, why is the default an alg that seems to give incorrect
> answers, instead of just defaulting to CDF?
> Also, is there a better way of doing this besides just appending CDF to the
> beginning every time?
>
> This is actually a real omission in the code. Using ComplextField(...)
instead of CDF is slower but provides more options, for instance working
with more than 53 bits. While computing eigenvalues of matrices numerically
is fundamentally a tricky operation that generally needs care and error
control, there are reasonable strategies to do a lot better than what
happens with ComplexField by default (as you can see with CDF). It would be
really good if matrices over ComplexField and RealField would use
algorithms that at least use some common sense. The warning about imprecise
answers would still apply, but to a much lower degree. (Technical note:
echelonization for CC and RR presently doesn't even use partial pivoting.
If we at least do that, we get something that's a little closer to the
"best effort" result people can reasonably expect.
Of course, by the time you're using numpy anyway, you should probably
prefer to use the state of the art numerical tools offered there. Only if
you need more than 53 bits should you consider the sage types.
--
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 post to this group, send email to sage-support@googlegroups.com.
Visit this group at https://groups.google.com/group/sage-support.
For more options, visit https://groups.google.com/d/optout.
