犰狳：arma::subview 的非连续视图

Question

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

即就像在

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

但是使用矩阵 A 的子视图？

即

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::subviewcs 作为函数参数的上下文。由于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));

因为子视图是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 或其他），但我宁愿不这样做，因为我已经使用犰狳多年而且我一直很满意...
• 我知道使用循环和更简单的子视图在我展示的更简单的代码中获得相同结果的可能性，但实际代码更复杂，并且此子视图将用于矩阵运算，所以我必须循环并复制临时对象中的子矩阵，我想避免这种情况

Answer 1

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

解决方案：没什么，我尝试使用 std::move() 将子视图的右值欺骗为空矩阵，对执行这 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) 是子视图的 r 值，而不是矩阵本身。