[Rcpp-devel] Detecting Singular Matrices in RcppArmadillo

[Alex Ilich]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20220128/05ebe0b5/attachment.html>