# Mathematica Precisions vs Doubles in C/C++

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

• "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. Apr 20, 2016 at 10:54
• But generally I would trust the values that agree with Mathematica's high precision result. Apr 20, 2016 at 10:56
• 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? Apr 20, 2016 at 11:10
• 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! Apr 20, 2016 at 12:06
• Can you give example matrices? Apr 20, 2016 at 12:55

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:= Eigenvalues[N[mat, 100]]

Out= {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:= Eigenvalues[N[mat]]

Out= {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:= << MATLink

In:= OpenMATLAB[]

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

Out= {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:= << LTemplate

SetDirectory[\$TemporaryDirectory];

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

code = "
#undef P

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:= CompileTemplate[tem,
"IncludeDirectories" -> {"/opt/local/include"},
"CompileOptions" -> {"-std=c++11", "-mmacosx-version-min=10.11"},
"LibraryDirectories" -> {"/opt/local/lib"},