Décomposition QR dans RcppArmadillo

Question

Vraiment confus pourquoi la sortie QR utilisant RcppArmadillo est différente de la sortie QR de R; La documentation Armadillo ne donne pas une réponse claire non plus. Essentiellement quand je donne à R une matrice Y qui est n * q (disons 1000 X 20), je récupère Q qui est 1000 X 20 et R 20 X 1000. C'est ce dont j'ai besoin. Mais quand j'utilise le solveur QR dans Armadillo, il me renvoie Q 1000 X 1000 et R 1000 X 20. Puis-je appeler la fonction qr de R à la place? J'ai besoin de Q pour avoir une dimension n x q, pas q x q. Le code ci-dessous est ce que j'utilise (c'est une partie d'une plus grande fonction).

Si quelqu'un peut suggérer comment le faire dans RcppEigen, ce serait utile aussi.

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

(NOTE: Cette réponse explique pourquoi les matrices R et retour RcppArmadillo avec des dimensions différentes, mais pas comment faire RcppArmadillo retourner la matrice 1000 * 20 que R ne Pour cela, peut-être jeter un oeil à la stratégie utilisée. près du fond de la définition de la fonction qr.Q().)

---

fonction de R qr() ne retourne pas directement Q. Pour cela, vous devez utiliser qr.Q(), comme ceci:

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 

Notez que qr.Q() retourne une matrice de la même dimension que m, au lieu d'un 5 complet * 5 matrice Q. Vous pouvez utiliser l'argument complete= pour contrôler ce comportement, et obtenir une matrice Q de dimension complète:

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 

Dans votre cas, il semble que RcppArmadillo est de retour la matrice complète 1000x1000 Q (comme qr.Q(q(m, complete=FALSE)) serait), plutôt que seulement ses 20 premières colonnes (comme qr.Q(q(m)) serait).