[Rcpp-devel] Detecting Singular Matrices in RcppArmadillo
Alex Ilich
alexilich93 at gmail.com
Sun Jan 30 17:11:31 CET 2022
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 <dmbates at 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 <alexilich93 at 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 <edd at 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 <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
>>>
>> _______________________________________________
>> 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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20220130/9b5e7e43/attachment.html>
More information about the Rcpp-devel
mailing list