使用 RcppArmadillo 修改输入

Question

我正在尝试使用 RcppArmadillo 通过完全旋转来实现 LU 分解。幸运的是，我有 this Matlab 代码可以满足我的要求，但我在将其转换为犰狳时遇到了一些挑战。

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

我知道使用常规 Rcpp 您可以像这样轻松修改输入：

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

这将 return 一个向量两倍的输入，并将输入更改为等于其原始值的两倍。但是，我很快发现这对 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 的 return 类型从 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)
*/

使用 RcppArmadillo 修改输入

Modifying inputs with RcppArmadillo

c++

r

armadillo

rcpp