我很困惑为什么使用RcppArmadillo输出的QR结果与R中的QR结果不同; Armadillo文档也没有给出清晰的答案。当我将一个n * q(假设1000 X 20)的矩阵Y输入R时,我得到了1000 X 20的Q和20 X 1000的R。这就是我需要的。但是当我在Armadillo中使用QR求解器时,它返回给我1000 X 1000的Q和1000 X 20的R。我能否调用R的qr函数?我需要Q的维度为n x q,而不是q x 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")
qr()
函数不能直接返回矩阵Q,需要使用qr.Q(qr(m))
来获取。返回的矩阵维度取决于complete=
参数的值。可以尝试以下代码来理解:m <- matrix(rnorm(10), ncol=2); qr.Q(qr(m)); qr.Q(qr(m), complete=TRUE)
。第一次调用qr.Q()
会返回一个5行2列的矩阵,而第二次则会返回完整的5行5列的Q矩阵,这取决于complete=
参数的设置。RcppArmadillo是否返回完整的1000行1000列的Q矩阵,而不是像qr.Q()
默认情况下只返回前20列的内容呢?(qr.R()
也有相同的complete=
参数)。 - Josh O'Brienqr.Q
中最后几行使用的策略。 - Josh O'Brien