使用犰狳访问下三角元素的更有效方法

Posted 2023-02-17

【中文标题】使用犰狳访问下三角元素的更有效方法

【英文标题】：More efficient way to access lower triangular elements using Armadillo

【发布时间】：2020-04-09 14:33:23

【问题描述】：

有没有更有效的方法来创建包含矩阵下三角元素的向量？对于某些算法，将这些元素放在一个向量中是很有用的。

但是，当我更喜欢指针时，下面的代码显然会创建副本。鉴于元素位置的非连续性质，这可能是不可能的。

我想到的另一种方法是通过find(trimatu(inmat)!=0) 左右创建一个索引矩阵，但我无法想象这样会更有效率。在双精度矩阵中找到精确的零通常不是很快，而且在我试图提取的三角形内也可能有实际的 0。 这里已经讨论了这些方面的内容（C++ Armadillo Access Triangular Matrix Elements），但是，这个问题已经存在 5 年了，从那时起犰狳已经有了很大的进步。

vec trimat2vec(mat const& inmat)

  int K = inmat.n_rows;
  int p = K*(K-1)/2;
  vec out(p);

  int counter=0;
  for(int i=0; i<K; i++)
    for(int j=1; j<K; j++)
      if(i<j)
        out(counter)=inmat(j,i);
        counter+=1;
      
    
  
  return(out);

【问题讨论】：

我不明白你为什么要使用向量。您将需要一个索引转换来从给定一些 (i,j) 索引的向量中获取正确的元素。您可以使用相同的转换（实际上是它的逆）来访问向量中的元素，因为它们在向量中的位置给定，而无需实际将元素复制到向量

只是为了方便。我可以对向量使用向量运算。我可以将矩阵中的所有信息存储在一个向量中，这很有帮助，因为犰狳只有立方体但没有 4d 数组。

好的，现在我明白了；）。我以为我可以在不了解犰狳的情况下提供帮助，但这似乎真的需要犰狳解决方案。 Fwiw，“稀疏矩阵”和“将矩阵展平为向量”将是我要寻找的。​​span>

是的，这是犰狳特有的。犰狳也很好地处理稀疏矩阵，但它甚至具有处理上/下三角矩阵的功能。从规则的方阵构造三角形很容易，但是，它不是从元素向量开始的。上面的代码sn-p也没有内置函数。

【参考方案1】：

好的，有了新版本的 Armadillo 9.870，现在有更好的方法来做到这一点。速度提升很小，尤其是在小矩阵上，但似乎是一致的。更重要的是，代码要短得多。

我调用新函数trimat2vec_new。

vec trimat2vec_new(mat const& inmat)
  return(inmat.elem(trimatl_ind(size(inmat),-1)));

R 中的简单基准测试：

a=matrix(1:900,30,30)

all(trimat2vec_new(a) == trimat2vec(a))

library(microbenchmark)
microbenchmark(trimat2vec_new(a),
               trimat2vec(a),
               times = 1000)

产量：

Unit: microseconds
              expr   min     lq     mean median    uq    max neval cld
 trimat2vec_new(a) 3.026 3.4060 3.633960  3.557 3.727 20.549  1000  a 
     trimat2vec(a) 3.116 3.6515 3.955116  3.787 3.958 42.981  1000   b

我特别喜欢新的trimatl_ind 功能的便利性。以及速度的小幅提升。