我想知道如何使用RcppArmadillo
. RcppArmadillo
使用以下方法在矩阵上运行 QR 分解相对简单:
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")
但是,当我将 的输出fun
与正常的 Rqr
方法进行比较时,我被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
我尝试在 R 源代码中查找旧的 Fortran,dqrdc2.f
但我不是一个强大的 Fortran 程序员,无法弄清楚这种紧凑形式是如何创建的。在这里和这里都提出了类似的问题,但我还没有找到解决方案。这里的最终目标是:
all.equal(fun(mat), qr(mat)$qr)
fun
如果我能了解这种紧凑的形式,显然会修改输出。
编辑
报告的文档qr
:
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.
这个“下三角”是我想了解的。