<div dir="ltr"><div>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.</div><div><br></div><div>Thanks,</div><div>Alex</div><div><br></div><div>library(Rcpp)<br>cppFunction(depends="RcppArmadillo", "NumericVector C_OLS(arma::mat X, arma::mat Y){<br>  arma::mat Xt = trans(X);<br>  arma::mat XtX = Xt * X;<br>  double d = det(XtX);<br>  Rcout << d; //print determinant<br>  if(d==0){<br>    NumericVector B2(X.n_cols, NA_REAL);<br>    return B2;<br>  } else{<br>    arma::mat XtX_inv= inv(XtX);<br>    NumericVector B = Rcpp::as<Rcpp::NumericVector>(wrap(XtX_inv * (Xt * Y)));<br>    return B;<br>  }}")<br><br>X<- matrix(data=c(-250, -250, -250, -250, -250, -250, 250, 200, 150,<br>                  100, 50, 0, 1, 1, 1, 1, 1, 1), ncol=3, nrow=6)<br><br><br>X2<- X<br><br>X2[1,2]<- X[1,2] + 5.684342e-14 #Add a slight numeric precision issue<br><br>Y<- matrix(1:6, ncol=1)<br><br>#Determinant of XtX in R is 0 for both<br>det(t(X) %*% X) ==0 #TRUE<br>det(t(X2) %*% X2) ==0 #TRUE<br><br>C_OLS(X,Y) #Determinant is 0 in C++ and correctly returns a vector of NA's<br>#0[1] NA NA NA<br><br>C_OLS(X2,Y) #Determinant is not 0 in C++ so the if statement is not triggered but the matrix is still uninvertible<br># 4.57764e-05<br># error: inv(): matrix is singular<br># Error in C_OLS(X2, Y) : inv(): matrix is singular</div></div>