[Rcpp-devel] Numerical precision in rotations with Eigen

[Dirk Eddelbuettel]

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.123234e-17)
| 
| I.e., the value for the Z component is 6.123234e-17. 
| 
| 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.776933e-13)
| 
| 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.220446e-16

>From the looks it seems like an instance of R FAQ 7.31:

   https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f

and the 'classic' referenced therein

   https://docs.oracle.com/cd/E19957-01/806-3568/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 1e-8.  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
| 
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel

-- 
dirk.eddelbuettel.com | @eddelbuettel | edd at debian.org