[Rcppdevel] Numerical precision in rotations with Eigen
Dirk Eddelbuettel
edd at debian.org
Wed Jan 10 13:22:37 CET 2024
On 10 January 2024 at 09:52, Rafael Ayala Hernandez wrote:
 Hi,

 I have implemented a function to rotate a 3D vector a given angle around a given axis (basically wrapping the functionality provided by Eigen::AngleAxis) as an Rcpp function.
 Below is an extract from the source file:

 #include <Rcpp.h>
 #include <RcppEigen.h>
 #include <Eigen/Eigen>
 #include <Eigen/Geometry>
 using namespace Rcpp;
 using Rcpp::as;
 using Eigen::Map;
 using Eigen::VectorXd;
 using Eigen::Vector3d;
 using Eigen::Matrix;
 using Eigen::MatrixXd;
 using Eigen::Matrix3d;
 using Eigen::Dynamic;
 using Eigen::AngleAxisd;

 // [[Rcpp::plugins("cpp11")]]
 // [[Rcpp::export]]
 NumericVector rotate3DVectorAngleAxis(NumericVector x, NumericVector axis, double angle) {
 // Note: x should be 1 single vector to rotate
 Map<VectorXd> xEigen(as<Map<VectorXd> >(x));
 Map<VectorXd> axisEigen(as<Map<VectorXd> >(axis));
 Matrix3d rotation;
 rotation = AngleAxisd(angle, axisEigen.normalized());
 Vector3d rotatedVector;
 rotatedVector = rotation.matrix() * xEigen;
 return wrap(rotatedVector);
 }


 The function executes with no problems and gives the expected output.

 However, I have noticed that there is times where a value of 0 would be expected for some vector components, but instead a very small value is returned. For example, executing the following on my machine:

 rotate3DVectorAngleAxis(c(0,0,1), c(1,0,0), pi/2)

 Which represents a rotation of 90 degrees around the X axis, and therefore the expected output value would be c(0,1,0)

 However, it instead results in an output vector of c(0.000000e+00, 1.000000e+00, 6.123234e17)

 I.e., the value for the Z component is 6.123234e17.

 I guess this is somehow related to machine precision, but is there an exact cause that could be possibly fixed?
 Additionally, why does the deviation from 0 only seem to happen for the Z component, and not for the X component?
 Varying the amount of rotation also seems to lead to a cumulative error on the same component. E.g., rotating by 4001 times pi/2 (i.e., 1000 complete 360 degrees rotations plus a 90 degree rotation):

 rotate3DVectorAngleAxis(c(0,0,1), c(1,0,0), 14401*pi/2)

 results in c(0.000000e+00, 1.000000e+00, 4.776933e13)

 Which seems again to be accumulating in the same Z component.

 Is there anything I could improve in my implementation to avoid these numerical errors?

 For reference, .Machine$double.eps on my machine is 2.220446e16
>From the looks it seems like an instance of R FAQ 7.31:
https://cran.rproject.org/doc/FAQ/RFAQ.html#Whydoesn_0027tRthinkthesenumbersareequal_003f
and the 'classic' referenced therein
https://docs.oracle.com/cd/E1995701/8063568/ncg_goldberg.html
It is probably your responsibility to identify such values and explicitly
round or set them to zero. I know nothing about the rotation here but in
other field (esp with iterations) a common rule of thumb is 'less than
sqrt(eps) is zero' meaning 1e8. Of course that is too general and may be
too broad as "it always depends".
Hope this helps a little. Maybe contact some field experts too.
Cheers, Dirk
 Thanks a lot in advance

 Best wishes,

 Rafa

 _______________________________________________
 Rcppdevel mailing list
 Rcppdevel at lists.rforge.rproject.org
 https://lists.rforge.rproject.org/cgibin/mailman/listinfo/rcppdevel

dirk.eddelbuettel.com  @eddelbuettel  edd at debian.org
More information about the Rcppdevel
mailing list