QR-декомпозиция в RcppArmadillo

Question

Действительно смущает, почему вывод QR с использованием RcppArmadillo отличается от QR-вывода от R; Документация Armadillo также не дает четкого ответа. По существу, когда я даю R матрицу Y, которая равна n * q (скажем, 1000 X 20), я возвращаюсь к Q, что составляет 1000 X 20 и R 20 X 1000. Это то, что мне нужно. Но когда я использую QR-решатель в Armadillo, он отбрасывает меня Q 1000 X 1000 и R 1000 X 20. Можно ли вместо этого вызвать функцию qr q? Мне нужно 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") 

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 возвращается полный 1000x1000 Q матрицы (как qr.Q(q(m, complete=FALSE)) бы), а не только его 1-й 20 столбцов (например, qr.Q(q(m))).