1
$\begingroup$

Normally, obtaining eigenvalues of random numerical matrices is fast. For instance a generic result looks like

Timing[Eigenvalues@RandomComplex[{-1 - I, 1 + I}, {500, 500}];]

{0.3125, Null}

Trying to evaluate the following takes so much time that I have to cancel the calculation

Timing[Eigenvalues@RandomComplex[{-1 - I, 1 + I}, {500, 500}, WorkingPrecision -> 16];]

The same is true even if I write:

Timing[Eigenvalues@RandomComplex[{-1 - I, 1 + I}, {500, 500}, WorkingPrecision -> $MachinePrecision];]

Why is this happening? And how can I still use Eigenvalues[] effectively when increasing precision?

$\endgroup$
  • 2
    $\begingroup$ You have forced a switch from machine to software arithmetic, and with it a change from calls to Lapack (using perhaps level 2 BLAS) to code that really cannot make fine-grained use of data locality or memory cache lines. $\endgroup$ – Daniel Lichtblau Aug 27 '15 at 14:44
  • $\begingroup$ @DanielLichtblau in "Some Notes on Internal Implementation" it says, "For dense arrays, LAPACK algorithms extended for arbitrary precision are used when appropriate." In the past I have tended to take this statement quite literally, but given your comment, it seems in doubt--or, at least, if LAPACK has been extended, BLAS has not been, even though it is a dependency of LAPACK. What is your point of view about this? $\endgroup$ – Oleksandr R. Aug 27 '15 at 15:27
  • 2
    $\begingroup$ @OleksandrR I believe that refers to Lapack code we rewrote for extended precision. The person who did that work is extremely well versed in Lapack intricacies, so it would be a faithful rewrite. Also it works at fixed precision so the tracking code will be largely disabled, hence should not be much of a factor. Or so I believe. $\endgroup$ – Daniel Lichtblau Aug 27 '15 at 16:43
  • $\begingroup$ @DanielLichtblau thanks. That is very useful. $\endgroup$ – Oleksandr R. Aug 27 '15 at 17:04
9
$\begingroup$

$MachinePrecision is different from MachinePrecision. The former calls for an arbitrary precision calcluation, done at the same precision as the machine-precision one. The main reason one would want to use this is to enable precision tracking, which is absent for a true machine-precision calculation using MachinePrecision.

And, there is your answer. Precision tracking, i.e. the augmentation of every basic arithmetic and higher-order mathematical operation to also determine whether it gains or loses precision, does not come for free. Twice as much information needs to be manipulated: both the actual numbers, and their remaining precision. And to calculate the change in precision, every operation implicitly has to be differentiated, so the number of basic arithmetic operations grows as well.

Functions that implement WorkingPrecision will tend also to increase the precision dynamically, if necessary, in order to ensure that the returned results have sufficient useful remaining precision. This is even more costly, and may be especially problematic in case any of your eigenvalues are very small, which could lead to large increases in the working precision in order to obtain any valid digits. You can disable this by setting $MinPrecision = $MaxPrecision = WorkingPrecision, but it will not stop precision from being tracked.

$\endgroup$
  • $\begingroup$ So there is no computationally efficient way to obtain all the eigenvalues at increased precision? It seems that the MachinePrecision number of digits is selected somewhat arbitrary to begin with... $\endgroup$ – Kagaratsch Aug 27 '15 at 16:32
  • 2
    $\begingroup$ Why do you say it is arbitrarily selected? It is based on the 53-bit significand of an IEEE754-compliant machine (64-bit) double precision quantity. I suppose it is arbitrary in the sense that we have 64-bit processors rather than, say, 71-bit ones, but there are reasons for that standardization. Anyway, I would think the answer is probably no, unless you can find some libraries that can calculate eigenvalues at arbitrary precision. It is not trivial to do, so I'm not sure how many such even exist outside of software like Mathematica and Maple. $\endgroup$ – Oleksandr R. Aug 27 '15 at 17:08
  • 2
    $\begingroup$ If 128-bit precision is good enough for you, you might find some libraries that support IEEE long doubles. $\endgroup$ – Oleksandr R. Aug 27 '15 at 17:12
  • $\begingroup$ Ok, I guess I should be looking for the 128-bit libraries then. $\endgroup$ – Kagaratsch Aug 27 '15 at 18:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.