[Rcpp-devel] Detecting Singular Matrices in RcppArmadillo

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

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

On 28 January 2022 at 13:25, Alex Ilich wrote:
| Thank you Neal and Dirk. Based on these comments I think using
arma::rcond

Ah, yes, of course. Should have remembered arma::rcond.

| to get the reciprocal of the condition factor, and checking if that is
zero
| rather than the determinant may be the route to go. At least for this
| example that works.
|
| 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 rcf = rcond(XtX);
|   Rcout << rcf; //print reciprocal of condition factor
|   if(rcf==0){

Re-read Neal's first email.  You probably do not want zero (apart from all
the R FAQ 7.31 reasons to not compare a double with equality :wink: )

Dirk

|     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)
|
| C_OLS(X,Y)
| # 0[1] NA NA NA
|
| C_OLS(X2,Y)
| # 0[1] NA NA NA
|
|
|
| On Fri, Jan 28, 2022 at 11:53 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
--
https://dirk.eddelbuettel.com | @eddelbuettel | edd at debian.org