QR decomposition in RcppArmadillo

Question

Really confused why the QR output using RcppArmadillo is different than QR output from R; Armadillo documentation doesnt give a clear answer either. Essentially when I give R a matrix Y that is n * q (say 1000 X 20 ) , I get back Q which is 1000 X 20 and R 20 X 1000. This is what I need. But when I use the QR solver in Armadillo, it throws me back Q 1000 X 1000 and R 1000 X 20. Can I call R's qr function instead? I need Q to have dimension n x q, not q x q. Code below is what I am using(its a part of a bigger function).

If someone can suggest how to do it in RcppEigen, that'd be helpful too.

library(inline)
library(RcppArmadillo)

src <- '
    Rcpp::NumericMatrix Xr(Xs);
    int q = Rcpp::as<int>(ys);

    int n = Xr.nrow(), k = Xr.ncol();
    arma::mat X(Xr.begin(), n, k, false);

    arma::mat G, Y, B;

    G = arma::randn(n,q);

    Y = X*G;

    arma::mat Q, R;
    arma::qr(Q,R,Y);

    return Rcpp::List::create(Rcpp::Named("Q")=Q,Rcpp::Named("R")=R,Rcpp::Named("Y")=Y);'


rsvd <- cxxfunction(signature(Xs="numeric", ys="integer"), body=src, plugin="RcppArmadillo")

R's `qr()` function does not return the `Q` matrix directly. Instead, you need to use `qr.Q(qr(m))`, and the dimension of the returned matrix will depend on the value of the `complete=` argument. Try this to see what I mean: `m <- matrix(rnorm(10), ncol=2); qr.Q(qr(m)); qr.Q(qr(m), complete=TRUE)`. The first call to `qr.Q()` returns a 5x2 matrix, while the second returns the full 5x5 Q matrix. Could it be that RcppArmadillo is returning the full 1000x1000 Q matrix, rather than just its 1st 20 columns (as `qr.Q()` would by default)? (`qr.R()` also has a `complete=` arg, for the same reason.) — Josh O'Brien, Jun 09 '12 at 19:57
You're right about using qr.Q(), I actually use it in the R-version. You're correct, RcppArmadillo returns the full 1000X1000. — pslice, Jun 09 '12 at 20:25
@JoshO'Brien +1 spot on. If you prepare that as an Answer that the OP can accept we can close this one off successfully. — Gavin Simpson, Jun 09 '12 at 20:28
I'm a little new to using RcppArmadillo/C++, so how I would I redefine Q after doing the QR such that it is just the first 20 columns? — pslice, Jun 09 '12 at 20:31
@pslice -- I haven't used RcppArmadillo or C++, so can't speak directly to that. If you want to do something other than make the full Q and then throw out all but the 1st 20 columns, maybe have a look at the strategy used in the last few lines of `qr.Q`. — Josh O'Brien, Jun 09 '12 at 20:35
@GavinSimpson -- OK, just did that, after rereading to see that I really had answered the main question. Thanks. — Josh O'Brien, Jun 09 '12 at 20:57

Josh O'Brien · Accepted Answer · 2012-06-09T21:18:06.727

(NOTE: This answer explains why R and RcppArmadillo return matrices with different dimensions, but not how to make RcppArmadillo return the 1000*20 matrix that R does. For that, perhaps take a look at the strategy used near the bottom of the qr.Q() function definition.)

R's qr() function does not return Q directly. For that, you need to use qr.Q(), like this:

m <- matrix(rnorm(10), ncol=2)
qr.Q(qr(m))
#             [,1]        [,2]
# [1,] -0.40909444  0.05243591
# [2,]  0.08334031 -0.07158896
# [3,]  0.38411959 -0.83459079
# [4,] -0.69953918 -0.53945738
# [5,] -0.43450340  0.06759767

Notice that qr.Q() returns a matrix of the the same dimension as m, rather than a full 5*5 Q matrix. You can use the complete= argument to control this behavior, and get a Q matrix of full dimension:

qr.Q(qr(m), complete=TRUE)
#             [,1]        [,2]       [,3]       [,4]        [,5]
# [1,] -0.40909444  0.05243591  0.3603937 -0.7158951 -0.43301590
# [2,]  0.08334031 -0.07158896 -0.8416121 -0.5231477  0.07703927
# [3,]  0.38411959 -0.83459079  0.2720003 -0.2389826  0.15752300
# [4,] -0.69953918 -0.53945738 -0.2552198  0.3453161 -0.18775072
# [5,] -0.43450340  0.06759767  0.1506125 -0.1935326  0.86400136

In your case, it sounds like RcppArmadillo is returning the full 1000x1000 Q matrix (like qr.Q(q(m, complete=FALSE)) would), rather than just its 1st 20 columns (like qr.Q(q(m)) would).

QR decomposition in RcppArmadillo

1 Answers1

Linked