3 added 4 characters in body
source | link

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).

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).

2 added 4 characters in body
source | link

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 CC++ code uses Armadillo.

I get near-identicalvery close but not identical results:

This seems to indicate that something might be wrong with your C code.

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 near-identical results:

This seems to indicate that something might be wrong with your C code.

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:

1
source | link

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 near-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}

This seems to indicate that something might be wrong with your C code.


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}