c++ - Armadillo:arma::subview 的非连续 View

Question

有没有办法使用矩阵的非连续 View

即喜欢在

arma::mat A = arma::zeros(3,3);
arma::uvec idx = {0,2};
A(idx,idx) += 2;

但使用矩阵 A 的 subview ？

即

arma::subview<double> swA = A.submat(0,0,2,2);
swA(idx,idx) += 2.5;

这最后一点不能编译为

 error: no match for call to ‘(arma::subview<double>) (arma::uvec&, arma::uvec&)’    
swA(idx,idx) += 2.5;

只是为了提供一些有助于我使用 arma::subviewc 作为函数参数的上下文。由于 A.submat(0,0,2,2) 是一个 rtype 我不能将它作为非常量参数传递给函数并且在函数内部我需要能够以非非连续的方式访问元素。

一个最小的(非工作的)例子也明白我的意思可能是下面的

#include <iostream>
#include <armadillo>

void f(arma::subview<double> x)
{ 
  arma::uvec i = {0,2};
  arma::uvec j = {1,2};
  x(i,j)  += 2.5;
}

int main ()
{

  arma::mat A = arma::zeros(5,5);
  std::cout << A << std::endl;

  f(A.submat(0,0,2,2));
  std::cout << A << std::endl;

  return 0;
}

gcc 再次返回的地方

 error: no match for call to ‘(arma::subview<double>) (arma::uvec&, arma::uvec&)’

解决这个问题的愚蠢方法是复制矩阵，将其作为引用传递，然后复制回 A:

#include <iostream>
#include <armadillo>

void f(arma::mat& x)
{ 
  arma::uvec i = {0,2};
  arma::uvec j = {1,2};
  x(i,j)  += 2.5;
}

int main ()
{

  arma::mat A = arma::zeros(5,5);
  std::cout << A << std::endl;

  arma::mat B = A.submat(0,0,2,2);
  f(B);
  A.submat(0,0,2,2) = B;
  std::cout << A << std::endl;

  return 0;
}

编译和运行完美，但我需要不惜一切代价避免复制矩阵(A 比 5x5 大得多!)

再次明确，我不能做以下事情

[...]
void f(arma::mat& x)
{ 
  arma::uvec i = {0,2};
  arma::uvec j = {1,2};
  x(i,j)  += 2.5;
}
[...]
f(A.submat(0,0,2,2));

因为 subview 是 rtype 并且我会从 gcc 获取

invalid initialization of non-const reference of type ‘arma::mat& {aka arma::Mat<double>&}’ from an rvalue of type ‘arma::mat {aka arma::Mat<double>}

我遇到麻烦了吗(只是一个实现问题，不在开发人员的 TODO 列表中)，还是比我聪明的人提供了一个优雅的解决方案？

谢谢!

旁注:

• 我愿意将线性代数库更改为其他库(例如 Eigen 或其他)，如果这样的事情在那里是微不足道的，但我宁愿不这样做，因为我多年来一直在使用 Armadillo ，而且我一直都很满意...
• 我知道可以使用循环和更简单的 subview 在我展示的更简单的代码中获得相同的结果，但实际代码更复杂，并且该 subview 将用于矩阵运算，所以我必须循环并复制临时对象中的子矩阵，我想避免这种情况

Answer 1

当您调用 subview 方法时，他/她/它使用不同的类 ( arma::matrix<T >operator() )在矩阵上....正如你发现的(和我一样，因为我来这里寻求答案)这是一个右值，所以它不能被引用......但是类的元素 arma::matrix<T>可以由类 subview 的对象初始化(那么为什么他们不创建一个访问函数来使用 subview 作为左值？为什么他们甚至创建类 subview ？)......我认为他们可以用一些模板和解决这个问题Move() , 但从底部继承到上层类，如 arma::mat . (你要求的问题不仅发生在稀疏矩阵上，甚至发生在密集矩阵上，只有你输入f(const arma:mat& )才有效，但你不能修改它)

解决方案:没什么，我尝试使用 std::move() 将 subview 的右值欺骗到一个空矩阵，对执行这 3 个循环 100000 次的时间进行基准测试；

/**********************************1************** ********/

void f(arma::mat &B)
{
}

// [[Rcpp::export]]
int func()
{
  arma::mat A = arma::zeros(50,50);
  //std::cout <<"A:\n"<< A << std::endl;

  arma::uvec i = arma::regspace<arma::uvec>( 0, 1, 10);
  arma::uvec j = arma::regspace<arma::uvec>( 0, 1, 10);
  unsigned p=1;
  do{
    arma::mat B=A(i,j);
    f(B);
    A(i,j)= B;
    p++;
  }while(p<100000);
  return 0;
}

/******************************2********************* ****/

void f(arma::mat &B)
{
}

// [[Rcpp::export]]
int func2()
{
  arma::mat A = arma::zeros(50,50);
  arma::uvec i = arma::regspace<arma::uvec>( 0, 1, 10);
  arma::uvec j = arma::regspace<arma::uvec>( 0, 1, 10);
  unsigned p=1;
  arma::mat B;
  do{
    B=A(i,j);
    f(B);
    A(i,j) = B;
    p++;
  }while(p<100000);
  return 0;
}

/***************************3************************ */

void f(arma::mat &B)
{
}

// [[Rcpp::export]]
int func3()
{
  arma::mat A = arma::zeros(50,50);
  arma::uvec i = arma::regspace<arma::uvec>( 0, 1, 10);
  arma::uvec j = arma::regspace<arma::uvec>( 0, 1, 10);
  unsigned p=1;
  arma::mat B;
  do{
    B=std::move(A(i,j));
    f(B);
    A(i,j)= std::move(B);
    p++;
  }while(p<100000);
  return 0;
}

从 R 编译和启动与构造函数复制或 std::move() 没有区别。甚至尝试使用 arma::mat B(std::move(A(i,j))没有区别，时间与 func() 相同(由于变量声明，略高于 func2() 和 func3())。

  Unit: milliseconds
    expr      min       lq     mean   median       uq      max neval cld
  func() 20.56159 20.68200 21.22679 20.98329 21.55746 27.62831  1000   b
 func2() 19.21450 19.37398 19.88782 19.66397 20.14731 27.86020  1000  a 
 func3() 19.19892 19.37547 19.89098 19.67835 20.17973 27.19525  1000  a 

这似乎是最快的(因为之前 A(i,j)= std::move(B) 将 B 大小设为 0 并且当 B= A(i,j) 时需要重新分配):

arma::mat B;
  do{
    B=std::move(A(i,j));
    f(B);
    A(i,j)= B;
    p++;
  }while(p<100000);
  return 0;

相同的基准

Unit: milliseconds
    expr      min       lq    mean   median       uq      max neval
 func4() 19.22223 19.30687 19.3907 19.33198 19.36293 23.20771  1000

std::move(A(i,j))似乎不影响(不松散的内容)矩阵 A，因为 A(i,j) 是 subview 的 r 值，而不是矩阵本身。