2

我想知道如何使用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.

这个“下三角”是我想了解的。

4

1 回答 1

1

我不明白 由 返回的紧凑形式qr,但您可以轻松地从对象中提取QR矩阵:

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

您会发现它们与 Rcpp 返回的内容相匹配。

我查看了手册,但我无法从中看出正面或反面。实际函数qr.Q只是调用一个 Fortran 函数,我无法解析。

于 2015-03-02T21:31:28.160 回答