在RcppArmadillo中的QR分解

Question

我很困惑为什么使用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")

- pslice

---

R的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'Brien

你关于使用qr.Q()是正确的，在R版本中我也用到了它。你说得对，RcppArmadillo返回完整的1000X1000矩阵。 - pslice

@JoshO'Brien +1，非常准确。如果您将其准备为答案，那么OP可以接受我们可以成功关闭此问题。 - Gavin Simpson

我对使用RcppArmadillo/C++还有点陌生，那么在进行QR分解后，如何重新定义Q使其仅包含前20列？ - pslice

@pslice -- 我没有使用过RcppArmadillo或C++，所以无法直接对此发表意见。如果您想要做的不仅仅是制作完整的Q矩阵，然后丢弃除了前20列之外的所有内容，也许可以看一下qr.Q中最后几行使用的策略。 - Josh O'Brien

@GavinSimpson -- 好的，我刚刚完成了这个任务，在重新阅读以确保我真正回答了主要问题后。谢谢。 - Josh O'Brien

Answer 1

注意：本答案解释了为什么R和RcppArmadillo返回的矩阵维度不同，但不解释如何使RcppArmadillo返回与R相同的1000*20矩阵。对于这个问题，请参考qr.Q()函数定义底部附近使用的策略。

---

R的qr()函数不直接返回Q。要获取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返回了完整的1000x1000 Q矩阵（就像qr.Q(q(m, complete=FALSE))一样），而不仅仅是它的前20列（就像qr.Q(q(m))一样）。

- Josh O'Brien