Разложение QR в RcppArmadillo

Действительно смущен, почему вывод 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. Могу ли я вместо этого вызвать функцию 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 rcpp decomposition armadillo
person pslice    schedule 09.06.2012    source источник
comment
Функция qr() языка R не возвращает напрямую матрицу 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, а второй возвращает полную Q-матрицу 5x5. Может быть, RcppArmadillo возвращает полную матрицу Q 1000x1000, а не только ее первые 20 столбцов (как qr.Q() по умолчанию)? (qr.R() также имеет аргумент complete= по той же причине.)   -  person Josh O'Brien    schedule 09.06.2012
comment
Вы правы насчет использования qr.Q(), я действительно использую его в R-версии. Вы правы, RcppArmadillo возвращает полные 1000X1000.   -  person pslice    schedule 10.06.2012
comment
@JoshO'Brien +1 место. Если вы подготовите это как ответ, который может принять ОП, мы сможем успешно закрыть его.   -  person Gavin Simpson    schedule 10.06.2012
comment
Я немного новичок в использовании RcppArmadillo/C++, поэтому как мне переопределить Q после выполнения QR, чтобы это были только первые 20 столбцов?   -  person pslice    schedule 10.06.2012
comment
@pslice - я не использовал RcppArmadillo или C++, поэтому не могу говорить об этом напрямую. Если вы хотите сделать что-то другое, а не составить полный Q, а затем выбросить все столбцы, кроме первых 20, возможно, взгляните на стратегию, использованную в последних нескольких строках qr.Q.   -  person Josh O'Brien    schedule 10.06.2012
comment
@GavinSimpson - Хорошо, только что сделал это, перечитав, чтобы убедиться, что я действительно ответил на главный вопрос. Спасибо.   -  person Josh O'Brien    schedule 10.06.2012

Ответы (1)


arrow_upward
6
arrow_downward

(ПРИМЕЧАНИЕ. Этот ответ объясняет, почему R и RcppArmadillo возвращают матрицы с разными размерами, но не объясняет, как заставить RcppArmadillo возвращать матрицу 1000 * 20, которую делает R. Для этого, возможно, взгляните на стратегию, используемую в нижней части функции qr.Q(). определение.)

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

person Josh O'Brien    schedule 09.06.2012