在RcppArmadillo QR分解(QR decomposition in RcppArmadi

2019-09-17 23:19发布

真糊涂为什么使用RcppArmadillo的QR输出大于来自R输出QR不同; 犰狳文档犯规要么给出一个明确的答案。 本质上讲,当我给R A矩阵Y是N * Q(比如1000×20),我得到的返回问题是1000×20和R 20×1000。这是我需要的。 但是,当我使用QR求解的犰狳,它抛出我返回问题1000×1000和R 1000×20.我能叫的r QR功能呢? 我需要Q可具有尺寸NXQ,不QX Q值。 下面的代码是什么我使用(它更大的功能的一部分)。

如果有人可以建议如何做到这一点的RcppEigen,那简直是太有帮助。

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")

Answer 1:

(注:此答案解释了为什么R和不同尺寸RcppArmadillo回报矩阵,而不是如何让RcppArmadillo返回1000 * 20矩阵的R确实对于这一点,也许看看附近的底部使用的策略。 qr.Q()函数的定义。)


的r qr()函数不直接返回Q值。 对于这一点,你需要使用qr.Q()就像这样:

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

注意, qr.Q()返回相同的尺寸的矩阵m ,而不是一个完整的5×5 Q矩阵。 您可以使用complete=参数来控制这种行为,并得到全尺寸的Q矩阵:

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

根据你的情况,这听起来像RcppArmadillo还是返回了1000×1000全Q矩阵(如qr.Q(q(m, complete=FALSE))会),而不仅仅是它的第一个20列(如qr.Q(q(m))会)。



文章来源: QR decomposition in RcppArmadillo