Thanks for the advice Douglas. Just using B = solve(X,Y) works and it outputs a warning that an approximate solution is attempted when the matrix is nearly singular. To suppress the warning I've just included #define ARMA_WARN_LEVEL 1 at the top of my cpp code. Also, no longer getting some odd outliers that I was previously! For documentation purposes, would it still be accurate to say "parameters are fit using ordinary least squares"?
Thanks, Alex On Sat, Jan 29, 2022 at 12:05 PM Douglas Bates <dmba...@gmail.com> wrote: > If you are concerned about singularity it would be better to use a QR or > singular value decomposition of X because the condition number of X'X is > the square of the condition number of X. Just the act of forming X'X from X > makes the conditioning much worse. > > I believe that arma::solve, as used by > https://github.com/RcppCore/RcppArmadillo/blob/master/src/fastLm.cpp, > uses a QR decomposition. > > In general it is much more stable to work with a decomposition of X than > to form X'X and work with it. > > Also, even if the formula for a calculation involves the inverse of a > matrix, it is rarely a good idea to evaluate the full inverse. Most of the > time, evaluating the inverse is much more work than is required to solve a > single linear system of equations. It is like saying that evaluating a > scalar quotient like a/b requires that you first evaluate the inverse of b, > 1/b, then form the product a * (1/b). In fact, it is worse because > evaluating the inverse of an n by n matrix A takes roughly the same amount > of work as solving n systems of equations of the form Ax = b. > > On Fri, Jan 28, 2022 at 6:56 PM <alexilic...@gmail.com> wrote: > >> 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 > >
_______________________________________________ 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