Mathematica Precisions vs Doubles in C/C++

Question

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

Answer 1

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}