CODEWITHSUNDEEP
c++precisioneigen3
Matt
Matt

Reputation: 2802

Numerical differences between Eigen version 3.2.0 and 3.3.4

I recently rolled a code that uses some basic vector math and SelfAdjointEigenSolver from Eigen 3.2.0 to 3.3.4. The code was prototyped with Matlab and it agreed with 3.2.0. The motivation for rolling to 3.3.4 was to avoid some warnings and errors that occurred with C++11(and newer) and 3.2.0.

The code works by reading thousands of points into a Eigen::Matrix<double, Eigen::Dynamic, 3>. When I inspect the values using gdb the input values are the same. One piece of code that shows a numerical differences is shown below. pilType is defined as double.

// this is the point cloud, center it about the origin
Eigen::Matrix<pilType, 1, 3> center = gridPnts.colwise().mean();
gridPnts.rowwise() -= center;
// Get the transformation matrix to align the point cloud with it's normal
// Build the co variance matrix
Eigen::Matrix<double, 3, 3> S = gridPnts.transpose() * gridPnts;
S /= static_cast<pilType>(gridPnts.rows() - 1);
Eigen::SelfAdjointEigenSolver<pilMat3> es(S);
Eigen::Matrix<pilType, 3, 2> trans;
trans = es.eigenvectors().block<3, 2>(0, 1);
// convert the point cloud to 2D for Convex Hull calculation
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> output(gridPnts.rows(), 2);
output = gridPnts * trans;

I am compiling the new version with GCC 5.3.0, flags -std=c++11 -O3 -DEIGEN_NO_DEBUG -Wall -Wextra -Werror -march=native -mtune=native. I have tried with and without the march and mtune options. The old version was compiled with GCC 4.3.4 and flags -O3 -DEIGEN_NO_DEBUG -Wall -Wextra -Werror. The CPU is a SandyBridge E5-2670 running SuSE Enterprise SLES 11. The differences start after the division of S by the number of points - 1. They are very small, on the order of 1e-12 but when it comes to the eigen solution the values very close to 0 can cause a change of sign in the eigenvectors. This causes the transformation matrix to be switched for one axis. The resulting output is flipped along this axis. The resulting 2D patch is sent to a routine to find the convex hull. The same points are returned but in reverse order due to the axis flip. This is not an error but it was unexpected.

In fact throughout the application various numbers show differences at the 1e-12 to 1e-15 range.

Here are some values when run with GDB (the top values are from version 3.3.4, the lower are from 3.2.0)

S before division:

    Eigen::Matrix<double,3,3,ColMajor> (data ptr: 0x7fffffffb0b0) = {[0,0] = 3532221.7020642869,
  [1,0] = 1.0913936421275139e-11, [2,0] = 3332.4628071428679, [0,1] = 1.0913936421275139e-11,
  [1,1] = 335265.83999999962, [2,1] = -1.3073986337985843e-12, [0,2] = 3332.4628071428679,
  [1,2] = -1.3073986337985843e-12, [2,2] = 4697.4509785714226}

Eigen::Matrix<double,3,3,ColMajor> (data ptr: 0x7fffffffb8f0) = {[0,0] = 3532221.7020642869,
  [1,0] = 1.0913936421275139e-11, [2,0] = 3332.4628071428679, [0,1] = 1.0913936421275139e-11,
  [1,1] = 335265.83999999962, [2,1] = -1.3073986337985843e-12, [0,2] = 3332.4628071428679,
  [1,2] = -1.3073986337985843e-12, [2,2] = 4697.4509785714226}

S after division

Eigen::Matrix<double,3,3,ColMajor> (data ptr: 0x7fffffffb0b0) = {[0,0] = 21151.028156073575,
  [1,0] = 6.5352912702246338e-14, [2,0] = 19.954867108639927, [0,1] = 6.5352912702246338e-14,
  [1,1] = 2007.5798802395186, [2,1] = -7.8287343341232597e-15, [0,2] = 19.954867108639927,
  [1,2] = -7.8287343341232597e-15, [2,2] = 28.128448973481571}

Eigen::Matrix<double,3,3,ColMajor> (data ptr: 0x7fffffffb8f0) = {[0,0] = 21151.028156073575,
  [1,0] = 6.5352912702246351e-14, [2,0] = 19.954867108639927, [0,1] = 6.5352912702246351e-14,
  [1,1] = 2007.5798802395188, [2,1] = -7.8287343341232597e-15, [0,2] = 19.954867108639927,
  [1,2] = -7.8287343341232597e-15, [2,2] = 28.128448973481575}

After preforming the eigensolution the value for trans are:

{[0,0] = 3.4096942364292876e-18, [1,0] = -1,
  [2,0] = 3.9893751393767394e-18, [0,1] = 0.99999955376919869, [1,1] = 3.4134614846087836e-18,
  [2,1] = 0.00094470175363752104}

{[0,0] = 0, [1,0] = 1,
  [2,0] = -3.2750362278258563e-15, [0,1] = 0.9999995537691988, [1,1] = 3.093932467653499e-18,
  [2,1] = 0.00094470175363752125}

The numerical differences are admittedly very small and are within the noise of accuracy. In the shown numbers the sign of the transformation matrix flipped which caused the 2D cloud to flip. Is there a way to clean up what was done such that these differences go away? I've seen some Eigen functions with thresholds for small values but I didn't see one for division.

Upvotes: 1

Views: 304

Answers (2)

Matt
Matt

Reputation: 2802

There was a change in the way that Eigen calculated division by scalar with version 3.2.7.

operator/=(Scalar) now performs a true division (instead of mat*(1/s))

Testing the code with S *= 1.0 / (static_cast<pilType>(gridPnts.rows() - 1)); shows that the numerical differences are resolved.

Upvotes: 0

ggael
ggael

Reputation: 29205

You should not rely on the sign on the eigenvectors at all (regardless of the solver and library), but rather adjust them to your needs. For instance, if you want to be sure to have a right-handed basis, then update the last one from the cross product of the first twos.

Upvotes: 2

Related Questions