RcppArmadillo中的QR分解答案

【问题标题】：QR decomposition in RcppArmadilloRcppArmadillo中的QR分解
【发布时间】：2012-06-13 09:50:48
【问题描述】：

真的很困惑为什么使用 RcppArmadillo 的 QR 输出与 R 的 QR 输出不同；犰狳文档也没有给出明确的答案。本质上，当我给 R 一个矩阵 Y 是 n * q (比如 1000 X 20 ) 时，我得到 Q 是 1000 X 20 和 R 20 X 1000。这就是我需要的。但是当我在犰狳中使用 QR 求解器时，它让我回到 Q 1000 X 1000 和 R 1000 X 20. 我可以调用 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")

【问题讨论】：

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() 返回一个 5x2 矩阵，而第二次调用返回完整的 5x5 Q 矩阵。会不会是 RcppArmadillo 返回完整的 1000x1000 Q 矩阵，而不仅仅是它的第 20 列（默认情况下 qr.Q() 会返回）？（qr.R() 也有一个 complete= arg，出于同样的原因。）
使用 qr.Q() 是对的，我实际上在 R 版本中使用它。你是对的，RcppArmadillo 返回完整的 1000X1000。
@JoshO'Brien +1 现场。如果您将其准备为 OP 可以接受的答案，我们可以成功关闭此答案。
我对使用 RcppArmadillo/C++ 有点陌生，那么在执行 QR 后我将如何重新定义 Q 使其只是前 20 列？
@pslice -- 我没有用过 RcppArmadillo 或 C++，所以不能直接说。如果你想做一些事情而不是做完整的 Q 然后扔掉除第一个 20 列之外的所有列，也许看看qr.Q 的最后几行中使用的策略。

标签： r rcpp decomposition armadillo

【解决方案1】：

（注意：这个答案解释了为什么 R 和 RcppArmadillo 返回具有不同维度的矩阵，而不是如何让 RcppArmadillo 返回 R 所做的 1000*20 矩阵。为此，也许可以看看在底部附近使用的策略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 正在返回完整的 1000x1000 Q 矩阵（如 qr.Q(q(m, complete=FALSE)) 会），而不仅仅是它的前 20 列（如 qr.Q(q(m)) 会）。

【讨论】：