Thanks! I'll use if(rcf <= std::numeric_limits<double>::epsilon()) instead of rcf==0.
On Fri, Jan 28, 2022 at 7:32 PM Dirk Eddelbuettel <e...@debian.org> wrote: > > 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 <alexilic...@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@lists.r-forge.r-project.org > | https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel > -- > https://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org >
_______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel