[Rcpp-devel] Detecting Singular Matrices in RcppArmadillo

[Neal Fultz]

There's an R object that has the machine precision, which could be a
reasonable threshold.

.Machine$double.eps

I believe there is a similar constant in the C++ standard library.

You might also try checking the condition of the matrix instead of the
determinant, but you might take a performance hit. EG:
https://www.quora.com/Why-is-a-condition-number-more-useful-than-a-determinant-for-spotting-problems-with-the-solution-to-ax-b

- Neal

On Fri, Jan 28, 2022 at 8:54 AM Alex Ilich <alexilich93 at gmail.com> wrote:

> Hi, I have some code I've written to calculate regression parameters using
> ordinarily least squares using RcppArmadillo. This is meant to be performed
> iteratively over many small subsets (sliding window over spatial raster
> data) so I'd like to automatically detect when the regression can't be
> performed and just return a vector of NA's. Currently I check to see if the
> determinant is 0 in order to see if the XtX matrix can be inverted and it
> seems to work in the vast majority of cases but I did run into a case where
> a small amount of precision error was introduced and for some reason this
> causes the determinant to be non-zero in C++ however, the matrix still is
> considered singular (uninvertible) leading to an error. I was wondering if
> anyone had any suggestions? Maybe there's a minimum value non-zero value
> where it's considered essentially singular that I could use instead of 0,
> or maybe there's a way to return a vector of NAs if an error is detected
> rather than stopping execution and printing an error message. I've also
> seen that pinv can be used, but from some tests I ran it's substantially
> slower, so I'd like to avoid using that. I've included some example code
> below. Any suggestions would be greatly appreciated.
>
> Thanks,
> Alex
>
> library(Rcpp)
> cppFunction(depends="RcppArmadillo", "NumericVector C_OLS(arma::mat X,
> arma::mat Y){
>   arma::mat Xt = trans(X);
>   arma::mat XtX = Xt * X;
>   double d = det(XtX);
>   Rcout << d; //print determinant
>   if(d==0){
>     NumericVector B2(X.n_cols, NA_REAL);
>     return B2;
>   } else{
>     arma::mat XtX_inv= inv(XtX);
>     NumericVector B = Rcpp::as<Rcpp::NumericVector>(wrap(XtX_inv * (Xt *
> Y)));
>     return B;
>   }}")
>
> X<- matrix(data=c(-250, -250, -250, -250, -250, -250, 250, 200, 150,
>                   100, 50, 0, 1, 1, 1, 1, 1, 1), ncol=3, nrow=6)
>
>
> X2<- X
>
> X2[1,2]<- X[1,2] + 5.684342e-14 #Add a slight numeric precision issue
>
> Y<- matrix(1:6, ncol=1)
>
> #Determinant of XtX in R is 0 for both
> det(t(X) %*% X) ==0 #TRUE
> det(t(X2) %*% X2) ==0 #TRUE
>
> C_OLS(X,Y) #Determinant is 0 in C++ and correctly returns a vector of NA's
> #0[1] NA NA NA
>
> C_OLS(X2,Y) #Determinant is not 0 in C++ so the if statement is not
> triggered but the matrix is still uninvertible
> # 4.57764e-05
> # error: inv(): matrix is singular
> # Error in C_OLS(X2, Y) : inv(): matrix is singular
> _______________________________________________
> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20220128/6a25fefe/attachment.html>