c++ - 使用 RcppArmadillo 修改输入

Question

我正在尝试使用 RcppArmadillo 通过完全旋转来实现 LU 分解。幸运的是我有 this可以执行我想要的操作的 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;
}

这将返回一个两倍于输入的 vector ，并将输入更改为等于其原始值的两倍。但是，我很快发现这对 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 个矩阵的最佳方法是什么？一份 list ？

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)
*/