<div dir="ltr">
<div>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"? <br></div><div><br></div><div>Thanks,</div><div>Alex</div>
</div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Sat, Jan 29, 2022 at 12:05 PM Douglas Bates <<a href="mailto:dmbates@gmail.com">dmbates@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>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.</div><div><br></div><div>I believe that arma::solve, as used by <a href="https://github.com/RcppCore/RcppArmadillo/blob/master/src/fastLm.cpp" target="_blank">https://github.com/RcppCore/RcppArmadillo/blob/master/src/fastLm.cpp</a>, uses a QR decomposition.<br></div><div><br></div><div>In general it is much more stable to work with a decomposition of X than to form X'X and work with it.</div><div><br></div><div> 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.</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, Jan 28, 2022 at 6:56 PM <<a href="mailto:alexilich93@gmail.com" target="_blank">alexilich93@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Thanks! I'll use if(rcf <= std::numeric_limits<double>::epsilon()) instead of rcf==0.</div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, Jan 28, 2022 at 7:32 PM Dirk Eddelbuettel <<a href="mailto:edd@debian.org" target="_blank">edd@debian.org</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br>
On 28 January 2022 at 13:25, Alex Ilich wrote:<br>
| Thank you Neal and Dirk. Based on these comments I think using arma::rcond<br>
<br>
Ah, yes, of course. Should have remembered arma::rcond.<br>
<br>
| to get the reciprocal of the condition factor, and checking if that is zero<br>
| rather than the determinant may be the route to go. At least for this<br>
| example that works.<br>
| <br>
| library(Rcpp)<br>
| <br>
| cppFunction(depends="RcppArmadillo", "NumericVector C_OLS(arma::mat X,<br>
| arma::mat Y){<br>
| arma::mat Xt = trans(X);<br>
| arma::mat XtX = Xt * X;<br>
| double rcf = rcond(XtX);<br>
| Rcout << rcf; //print reciprocal of condition factor<br>
| if(rcf==0){<br>
<br>
Re-read Neal's first email. You probably do not want zero (apart from all<br>
the R FAQ 7.31 reasons to not compare a double with equality :wink: )<br>
<br>
Dirk<br>
<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 *<br>
| 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>
| C_OLS(X,Y)<br>
| # 0[1] NA NA NA<br>
| <br>
| C_OLS(X2,Y)<br>
| # 0[1] NA NA NA<br>
| <br>
| <br>
| <br>
| On Fri, Jan 28, 2022 at 11:53 AM Alex Ilich <<a href="mailto:alexilich93@gmail.com" target="_blank">alexilich93@gmail.com</a>> wrote:<br>
| <br>
| > Hi, I have some code I've written to calculate regression parameters using<br>
| > ordinarily least squares using RcppArmadillo. This is meant to be performed<br>
| > iteratively over many small subsets (sliding window over spatial raster<br>
| > data) so I'd like to automatically detect when the regression can't be<br>
| > performed and just return a vector of NA's. Currently I check to see if the<br>
| > determinant is 0 in order to see if the XtX matrix can be inverted and it<br>
| > seems to work in the vast majority of cases but I did run into a case where<br>
| > a small amount of precision error was introduced and for some reason this<br>
| > causes the determinant to be non-zero in C++ however, the matrix still is<br>
| > considered singular (uninvertible) leading to an error. I was wondering if<br>
| > anyone had any suggestions? Maybe there's a minimum value non-zero value<br>
| > where it's considered essentially singular that I could use instead of 0,<br>
| > or maybe there's a way to return a vector of NAs if an error is detected<br>
| > rather than stopping execution and printing an error message. I've also<br>
| > seen that pinv can be used, but from some tests I ran it's substantially<br>
| > slower, so I'd like to avoid using that. I've included some example code<br>
| > below. Any suggestions would be greatly appreciated.<br>
| ><br>
| > Thanks,<br>
| > Alex<br>
| ><br>
| > library(Rcpp)<br>
| > cppFunction(depends="RcppArmadillo", "NumericVector C_OLS(arma::mat X,<br>
| > 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 *<br>
| > 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<br>
| > 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<br>
| ><br>
| _______________________________________________<br>
| Rcpp-devel mailing list<br>
| <a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank">Rcpp-devel@lists.r-forge.r-project.org</a><br>
| <a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" rel="noreferrer" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
-- <br>
<a href="https://dirk.eddelbuettel.com" rel="noreferrer" target="_blank">https://dirk.eddelbuettel.com</a> | @eddelbuettel | <a href="mailto:edd@debian.org" target="_blank">edd@debian.org</a><br>
</blockquote></div>
_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" rel="noreferrer" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a></blockquote></div>
</blockquote></div>