Eigen floating-point precision

Question

I am having trouble with floating-point precision by using Eigen.

I have two Eigen::MatrixXd; the first matrix A (nx1) contains only positive integers numbers, while the second matrix B (nx1) contains a single column filled with the same real number (ex: -0.714312).

I need to compute the following Eigen::MatrixXd:

const auto exponential = [](double x)
{ return std::exp(x); };


MatrixXd W = B.unaryExpr(exponential);
MatrixXd residuals = A - W;

The problem is that when I print the sum of the residuals:

cout << residuals.sum();
// output = 6.16951e-06

I get a different value that by performing the same operation using R and the same input matrices.

By using R matrices I get -2.950208e-09. While the sum of the elements of A, B and W are the same both in C++ and in R.

Retch

Answer 1

It could be that R is using the x87 FPU with extended precision (80bits) whereas Eigen is using SSE units (64bits/double). You can check that by using Matrix<long double,Dynamic,Dynamic> matrix types or by making sure that your compiler will target the x87 FPU units.

Goaltender

Unfortunately using long double gives the same result. How can I check that my compiler targets the x87 FPU units? – Retch

Then please be more specific with actual values for A and B so that we can reproduce. – Goaltender