2

I was wondering how to return the 'compact' form from QR decomposition using RcppArmadillo. It is relatively simple to run a QR decomposition on a matrix using RcppArmadillo with the following:

library(inline)

src <- '
using namespace Rcpp;
using namespace arma;
mat X = as<mat>(A);
mat Q, R;

qr(Q,R,X);

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

fun <- cxxfunction(signature(A = "matrix"), src, plugin="RcppArmadillo")

However, when I compare the output of fun to the normal R qr method I am stumped by the 'compact' matrix returned by qr$qr.

mat <- matrix(seq(9), 3)

# Armadillo method
fun(mat)
$Q
           [,1]       [,2]       [,3]
[1,] -0.2672612  0.8728716  0.4082483
[2,] -0.5345225  0.2182179 -0.8164966
[3,] -0.8017837 -0.4364358  0.4082483

$R
          [,1]      [,2]          [,3]
[1,] -3.741657 -8.552360 -1.336306e+01
[2,]  0.000000  1.963961  3.927922e+00
[3,]  0.000000  0.000000  2.220446e-15

# normal R method
qr(mat)$qr
           [,1]      [,2]          [,3]
[1,] -3.7416574 -8.552360 -1.336306e+01
[2,]  0.5345225  1.963961  3.927922e+00
[3,]  0.8017837  0.988693  1.776357e-15

I have tried looking a the old Fortran in the R source dqrdc2.f but I am not a strong Fortran programmer and cannot figure out how this compact form is created. Similar questions were asked here and here but I haven't found a solution to this. The ultimate goal here is to have:

all.equal(fun(mat), qr(mat)$qr)

where clearly the output of fun would be modified if I can learn about this compact form.

EDIT

The documentation for qr reports:

Value

The QR decomposition of the matrix as computed by 
LINPACK or LAPACK. The components in the returned value correspond 
directly to the values returned by DQRDC/DGEQP3/ZGEQP3.

qr  
a matrix with the same dimensions as x. The upper triangle contains the 
\bold{R} of the decomposition and the lower triangle contains information 
on the \bold{Q} of the decomposition (stored in compact form). Note that 
the storage used by DQRDC and DGEQP3 differs.

This "lower triangle" is what I want to understand.

Community
  • 1
  • 1
cdeterman
  • 19,630
  • 7
  • 76
  • 100

1 Answers1

1

I don't understand the compact form returned by qr, but you can easily just extract the Q and R matrices from the object:

qr.Q(qr(mat))
qr.R(qr(mat))

You'll find that they match what is returned by Rcpp.

I looked at the manual, but I can't make heads or tails out of it. The actual function for qr.Q just calls a Fortran function, which I'm not able to parse.

nograpes
  • 18,623
  • 1
  • 44
  • 67
  • Thank you for your response, I was aware of those two functions but I was curious how the compact form was made. I will leave this question open to see if anyone does know. – cdeterman Mar 02 '15 at 21:33