用RcppArmadillo修改输入

Question

我正在尝试使用RcppArmadillo实现完全旋转的LU分解。幸运的是，我有这 Matlab代码来做我想做的事情，但是我在将它转换成Armadillo时遇到了一些困难。

我想让我的gecpLU函数像arma::LU一样工作，输入L、U和P，arma::LU函数修改输入的L、U和P矩阵，而不是返回L、U和P。

我知道，使用常规的Rcpp，您可以很容易地修改输入，如下所示：

NumericVector example(NumericVector X) {
    X = 2 * X;
    return X;
}

这将返回一个向量两倍于输入，并将输入更改为等于其原始值的两倍。但是，我很快发现这对RcppArmadillo不起作用。

arma::colvec example(arma::colvec X) {
    X = 2 * X;
    return X;
}

我知道当暴露在R中时，这不会改变输入，因为arma对象是R对象的副本，所以不能直接修改R对象，但我觉得我仍然应该能够编写类似Armadillo的LU这样的函数：

#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppArmadillo)]]
int gecpLU(arma::mat L, arma::mat U, arma::mat P, arma::mat Q, arma::mat A) {
  // Take A and overwrite LUPQ such that A=P*L*U*Q
  int n=A.n_rows;
  P.eye(n,n);
  Q.eye(n,n);
  arma::mat AA=A;
  // for (int i=0;i<(n-1);i++) {
  // delete a whole bunch of stuff not relevant to question
  // }
  L.eye(n,n);
  arma::mat tempmat=arma::trimatl(AA);
  tempmat.diag()*=0;
  L=L-tempmat;
  U=arma::trimatu(AA);

  return 0;
}

// [[Rcpp::export]]
List test(arma::mat A) {
  arma::mat L1,U1,P1,L2,U2,P2,Q;
  arma::lu(L1,U1,P1,A);
  gecpLU(L2,U2,P2,Q,A);
  return List::create(_["L1"]=L1,
                      _["U1"]=U1,
                      _["P1"]=P1,
                      _["L2"]=L2,
                      _["U2"]=U2,
                      _["P2"]=P2,
                      _["Q"]=Q);
}

在这种情况下，我不是将R矩阵传递给我的gecpLU函数，而是arma::mat，因此它应该能够修改输入。

当我运行test时，我得到了L1、U1和P1的矩阵，而L2、U2、P2和Q的0x0矩阵，我觉得我一定是误解了什么。可以用RcppArmadillo修改输入吗？如果不是，输出4个矩阵的最佳方法是什么？一份名单？

Answer 1

下：

int gecpLU(arma::mat L, arma::mat U, arma::mat P, arma::mat Q, arma::mat A)

您正在创建每个矩阵的新副本，然后在函数的末尾销毁它们。

您期望的是对象会被修改。为此，需要在对象类型的末尾追加一个&，以便编译器知道如何获得对象的引用。

void gecpLU(arma::mat& L, arma::mat& U, arma::mat& P, arma::mat& Q, arma::mat& A)

注意，我还将gecpLU的返回类型从int更改为void。请参见：

#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppArmadillo)]]
void gecpLU(arma::mat& L, arma::mat& U, arma::mat& P, arma::mat& Q, arma::mat& A) {
  // Take A and overwrite LUPQ such that A=P*L*U*Q
  int n=A.n_rows;
  P.eye(n,n);
  Q.eye(n,n);
  arma::mat AA=A;
  // for (int i=0;i<(n-1);i++) {
  // delete a whole bunch of stuff not relevant to question
  // }
  L.eye(n,n);
  arma::mat tempmat=arma::trimatl(AA);
  tempmat.diag()*=0;
  L=L-tempmat;
  U=arma::trimatu(AA);
}

// [[Rcpp::export]]
List test(arma::mat A) {
  arma::mat L1,U1,P1,L2,U2,P2,Q;
  arma::lu(L1,U1,P1,A);
  gecpLU(L2,U2,P2,Q,A);
  return List::create(_["L1"]=L1,
                      _["U1"]=U1,
                      _["P1"]=P1,
                      _["L2"]=L2,
                      _["U2"]=U2,
                      _["P2"]=P2,
                      _["Q"]=Q);
}

简单的临时示例，演示引用传递(允许修改)与按副本传递(最终破坏对象)。

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;

void reference_obj(arma::vec& y){
  y.fill(1);
}

void copy_obj(arma::vec y){
  y.fill(0);
}

// [[Rcpp::export]]
arma::vec on_reference_mod(arma::vec x) {

  reference_obj(x);

  return x;
}


// [[Rcpp::export]]
arma::vec on_copy_mod(arma::vec x) {

  copy_obj(x);

  return x;
}

/*** R
# Should get a vector of 1's
on_reference_mod(1:10)

# Should get a vector of 1:10
on_copy_mod(1:10)
*/