3
$\begingroup$

I'm having a bit of an issue regarding numerical precision and I'm not sure how to deal with it.

I have a certain randomly generated matrix, say $M$, that I wish to compute the eigenvalues. The entries of this matrix are real and span many orders of magnitude, from $10^{16}$ down to $10^{-14}$.

I am particularly interested in the case where the smallest eigenvalue is bounded in a region, say $10^{-10}<\lambda_{min}<10^{-9}$. My program then generates some random numbers for the entries of $M$, computes the eigenvalues through

EigenValues[N[M,precision]]

where 'precision' is specified somewhere before the above evaluation$^a$, and then stores some of the parameters for $M$ that meets the criteria above.

In order to test the numerical precision of my numbers I generated an example instance of $M$, let's call it $M0$. I write down its entries in exact form, i.e.

Precision[M0]=Infinity

I then compute the eigenvalues with different precisions, namely consider the two cases

Eigenvalues[N[M0,MachinePrecision]]
Eigenvalues[N[M0,30]]

and test these results against the same computation carried out by a C program.

For a specific matrix $M0$ that I am using for this test I get the following results for its eigenvalues.

  • Mathematica With MachinePrecision: $$ 1.08568\times 10^{15} \\ -1.08568\times 10^{15} \\ 2.04979\times 10^{13} \\ 1.98037\times 10^{13} \\ -1.51195\times 10^{10} \\ 1.51195\times 10^{10} \\ 737.092 \\ 709.496 \\ -288.09 \\ 287.764 \\ -0.0527701 $$

  • Mathematica With precision=30: $$ 1.08568468076658564504709826065\times 10^{15} \\ -1.08568468076544177351940777511\times 10^{15} \\ 2.04978825179076415641037942687\times 10^{13} \\ 1.98037159026622180001959563609\times 10^{13} \\ -1.51195073626019461804045801802\times 10^{10} \\ 1.51195073626019461804034605668\times 10^{10} \\ 737.067283246721413545966471801 \\ 709.473111591071399129829176386 \\ -288.009491702297922120741501773 \\ 287.845096874967906231858249555 \\ -\text{1.04616771734909972710824504576458020590918184213753$\grave{ }$30.*${}^{\wedge}$-8} $$

  • Custom made C code, where the numbers are taken to be Double $$ +2.472\times10^{-10}\\ +2.878\times10^{+02}\\ -2.880\times10^{+02}\\ +7.095\times10^{+02}\\ +7.371\times10^{+02}\\ +1.512\times10^{+10}\\ -1.512\times10^{+10}\\ +1.980\times10^{+13}\\ +2.050\times10^{+13}\\ -1.086\times10^{+15}\\ +1.086\times10^{+15} $$

Now, reading through Mathematica's documentation, I was expecting the custom made code and MachinePrecision to agree since

The typical arrangement is that all machine‐precision numbers in the Wolfram Language are represented as "double‐precision floating‐point numbers" in the underlying computer system. On most current computers, such numbers contain a total of 64 binary bits, typically yielding 16 decimal digits of mantissa.

And I understand that due to a mantissa of around 16 digits, the first numerical problems would start at around $0.01$, since the biggest number is of order $10^{15}$. This expectation is confirmed by the differences between the numbers obtained within Mathematica with different precisions. But I cannot make sense of differences between the C code and MachinePrecision numbers obtained in Mathematica.

My questions are then:

What numbers should I trust?

Why is MachinePrecision result different than C-code result?

Notes:

a) 'precision' is used as a defining variable whenever I can set Precision or WorkingPrecision. This includes the RandomReal[] functions calls, N[], etc. A numerical generated matrix $M$ has

Precision[M]=precision
$\endgroup$
10
  • 8
    $\begingroup$ "and test these results against the same computation carried out by a C program" <- There are different ways to compute eigenvalues, and they may give different results in precision-sensitive situations. What does your C program do exactly? Mathematica can also use different algorithms and unfortunately not all of these are fully documented. $\endgroup$
    – Szabolcs
    Apr 20, 2016 at 10:54
  • $\begingroup$ But generally I would trust the values that agree with Mathematica's high precision result. $\endgroup$
    – Szabolcs
    Apr 20, 2016 at 10:56
  • 6
    $\begingroup$ You didn't mention what exactly is done in your "custom-made C code", so it's hard to say anything. For instance, unsymmetric matrices whose entries span many orders of magnitude are usually preprocessed by a balancing transformation before being fed to a QR routine, to use a typical workflow example. Did your code do that? $\endgroup$ Apr 20, 2016 at 11:10
  • $\begingroup$ Hi guys, thanks for your comments! I'm working with other people and I didn't write the C code, but as far as I know it is a C++ code that uses Armadillo package. The matrix is Symmetric, not too big (11x11), with doubles as entries. And I assume it was compiled in a 64bit CPU. The fact that it is Symmetric from the start should make the usage of eigenvalue function straightforward, no? It has been a long time since my Numerical Analysis university course! $\endgroup$
    – romanovzky
    Apr 20, 2016 at 12:06
  • $\begingroup$ Can you give example matrices? $\endgroup$
    – Szabolcs
    Apr 20, 2016 at 12:55

1 Answer 1

3
$\begingroup$

This is just a long comment really.

Out of curiosity, I compared the result of three packages: Mathematica, MATLAB and eig_sym from Armadillo (compiled on OS X). You said that your C++ code uses Armadillo.

I get very close but not identical results:

Mathematica MachinePrecision:

{1.08568*10^15, -1.08568*10^15, 2.04979*10^13, 
 1.98037*10^13, -1.51195*10^10, 
 1.51195*10^10, 737.092, 709.496, -288.09, 287.764, -0.0531685}

MATLAB:

{1.08568*10^15, -1.08568*10^15, 2.04979*10^13, 
 1.98037*10^13, -1.51195*10^10, 
 1.51195*10^10, 737.092, 709.496, -288.09, 287.764, -0.0531657}

Armadillo eig_sym:

{1.08568*10^15, -1.08568*10^15, 2.04979*10^13, 
 1.98037*10^13, -1.51195*10^10, 
 1.51195*10^10, 737.092, 709.496, -288.089, 287.765, -0.0532012}

These are all different form the result obtained with Mathematica arbitrary precision arithmetic with 100 digits, which gives -1.046*10^-8 for the smallest eigenvalue.

Your main question is: which result can be trusted. We can try to compute the eigenvectors as well and verify:

{val, vec} = Eigensystem[N[mat, 100]];

mat.vec[[-1]]/vec[[-1]]

This gives a vector where all entries agree to many digits and they are all -1.046*10^-8.

Now let us try machine precision for the eigensystem computation, but use high precision again for the verification:

{val, vec} = Eigensystem[N[mat, MachinePrecision]];

v = SetPrecision[vec[[-1]], 100];

mat.v/v

{742.726, 704.17, 
 3.00435*10^8, 48.2152, -129.252, 39.4756, -6.75853*10^-8, \
-2.58842*10^6, 1.96944*10^13, 1.98037*10^13, -1.67772*10^7}

Now we don't get the same values.

This indicates that the machine precision result is incorrect, while the high precision result is correct (as one would naively expect).


The code I used, for reference:

In[2]:= Eigenvalues[N[mat, 100]]

Out[2]= {1.\
0856846807665856450470982606468499998443593068923500101675873225220935\
64695808687498260652176359947*10^15, \
-1.0856846807654417735194077751054318897882403216356104198311574238345\
32051036384281450170890942140023*10^15, 
 2.0497882517907641564103794268674566843581943794938406285013154951144\
39397586760816469279826548743094*10^13, 
 1.9803715902662218000195956360852879283233255266649313821701000225618\
94167023214407217436135853211420*10^13, \
-1.5119507362601946180404580180226453265713175363317940722790080527023\
38252239086243921340499166490296*10^10, 
 1.5119507362601946180403460566794916459982630836031304301486393770875\
41052713765641281269928732700045*10^10, \
737.067283246721413545966471801175928892427820125177205599257512879741\
5217269062183664175848994015722, \
709.473111591071399129829176386128726837373250051028953864335805315683\
3282193626859167101643402616852, \
-288.00949170229792212074150177336349486164509410462852414839291990667\
67324980503281076889597723655422, \
287.845096874967906231858249554911428826683897345958518850569903869812\
2959065307646083681528989497732, \
-1.0461677173490997271082450457645821365852016628481886960614069559648\
94667505277812120367633524821533*10^-8}

In[3]:= Eigenvalues[N[mat]]

Out[3]= {1.08568*10^15, -1.08568*10^15, 2.04979*10^13, 
 1.98037*10^13, -1.51195*10^10, 
 1.51195*10^10, 737.092, 709.496, -288.09, 287.764, -0.0531685}

In[4]:= AddPath["MATLink1Dev"]

In[7]:= << MATLink`

In[8]:= OpenMATLAB[]

In[11]:= Reverse@SortBy[Abs]@Flatten[MFunction["eig"][mat]]

Out[11]= {1.08568*10^15, -1.08568*10^15, 2.04979*10^13, 
 1.98037*10^13, -1.51195*10^10, 
 1.51195*10^10, 737.092, 709.496, -288.09, 287.764, -0.0531657}

In[12]:= << LTemplate`

SetDirectory[$TemporaryDirectory];

In[14]:= tem = LClass["Eig",
   {LFun["eig", {{Real, 2, "Constant"}}, {Real, 1}]}
   ];

code = "
  #undef P

  #include <armadillo>

  struct Eig {
    mma::RealTensorRef eig(mma::RealMatrixRef t) {
        arma::mat m( t.begin(), t.cols(), t.rows() );
        arma::vec ev = eig_sym(m.t());
        return mma::makeVector<double>(ev.size(), ev.memptr());
    }
  };
  ";
Export["Eig.h", code, "String"];

In[40]:= CompileTemplate[tem, 
 "IncludeDirectories" -> {"/opt/local/include"}, 
 "CompileOptions" -> {"-std=c++11", "-mmacosx-version-min=10.11"}, 
 "LibraryDirectories" -> {"/opt/local/lib"}, 
 "Libraries" -> {"armadillo"}]

In[41]:= LoadTemplate[tem]

eig = Make[Eig];

In[45]:= Reverse@SortBy[Abs]@eig@"eig"[mat]

Out[45]= {1.08568*10^15, -1.08568*10^15, 2.04979*10^13, 
 1.98037*10^13, -1.51195*10^10, 
 1.51195*10^10, 737.092, 709.496, -288.089, 287.765, -0.0532012}
$\endgroup$
3
  • $\begingroup$ Thanks a lot for this answer. I just want to cross-check the C++ results myself first before I give you a tick for correct answer :) $\endgroup$
    – romanovzky
    Apr 22, 2016 at 13:37
  • $\begingroup$ I got quite different results, mainly the lowest eigenvalue I got 1.0785e-06, using eig_sym. Using eig_gen I get -4.130e-08 which again is not compatible. I then decided to implement the check m.vec^T-vec^T.diag(vals) and could see that C and MachinePrecision results are not reliable, while setting precision to rel. high values provides consistent results. Long story short, I'll keep to Mathematica as I can be sure the precision is being tracked correctly. Thank you a lot! $\endgroup$
    – romanovzky
    Apr 25, 2016 at 10:16
  • $\begingroup$ @romanovzky I don't know much about linear algebra software, but I think that Armadillo's results may be influenced by the underlying BLAS/LAPACK implementations. I'm on OS X using it's built-in Accelerate framework (the default for Armadillo). On Linux Armadillo would be using some other BLAS depending on its configuration. $\endgroup$
    – Szabolcs
    Apr 25, 2016 at 10:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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