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.