非连续矩阵的高级构造函数

Posted 2023-02-16

【中文标题】非连续矩阵的高级构造函数

【英文标题】：Advanced constructor for non-contiguous matrices

【发布时间】：2019-05-31 05:22:04

【问题描述】：

在我的实现中，我经常使用子矩阵和矩阵块。我想知道犰狳中是否有一种方法可以让我提取一个更大的矩阵块，并为这个子矩阵使用与原始矩阵中的块相同的内存。我的问题是我不知道该怎么做，因为原始矩阵中的位置不连续。

这是一个简单的例子，说明当我的原始矩阵是 A = [A1 A2] 时我想要做什么：

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat foo(arma::mat A, arma::uword block_start, arma::uword block_length) 
  arma::uword rows = A.n_rows;
  arma::mat B = arma::mat(A.begin_col(block_start), rows, block_length, false);
// fill B here for illustration; 
// in the real piece of code I do various manipulations, multiplications, etc using B
  B.fill(1.0);
  return A;


/*** R
A <- matrix(0, 4, 4)
foo(A, 0, 2)
> A <- matrix(0, 4, 4)
> foo(A, 0, 2)
     [,1] [,2] [,3] [,4]
[1,]    1    1    0    0
[2,]    1    1    0    0
[3,]    1    1    0    0
[4,]    1    1    0    0
*/

在这种情况下，子矩阵的位置是连续的，我可以使用高级构造函数来链接内存。

现在假设我希望子矩阵为A[1:2, 1:2]。我可以在犰狳中获得一份B 的副本，它使用与A 中的原始元素相同的内存吗？ （理想情况下，这个问题的解决方案也可以推广到列也不连续的情况，例如A[c(1, 3), c(1, 3)]。）

编辑：澄清一下，我真的需要上面函数中的矩阵B 独立存在。我不在我的真实代码中使用fill 它，而是在各种矩阵乘法等中使用它（和多个其他子矩阵）。所以我追求的是一种将B 创建为@ 的非连续子矩阵的方法987654330@ 同时确保它们使用相同的内存。

【问题讨论】：

我想你想要一个圆孔中的方形钉。想想 5x5 矩阵。现在您想要内部 3x3 ... 但如果没有副本就不能连续拥有它，因为外部矩阵的大小决定了 5x5 中 25 个单元的使用方式。

@DirkEddelbuettel 你可能是对的。但这都是关于指针的，不是吗？一切都可以指向。理想情况下，可以使用submat(rows, cols).begin()，但这不起作用。当然，这可能会更不安全，因为你有更复杂的依赖。我觉得一次又一次地提取和复制子矩阵有点脏，所以避免这种情况会很棒。

嗯，我担心你只是在重复你想要的东西，而忽略了内存的使用方式。从记事本开始，“绘制”一个 5x5 矩阵。然后画出第一个、第二个、..、第 20 个地址是什么。然后尝试找到一个多行/多列子视图并说服自己它不能是连续的子向量。我可能是错的，但我仍然认为你在要求一个不存在的东西。

@DirkEddelbuettel 我明白你的意思。仅仅证明我所要求的不可行对我来说很有价值。缺少功能通常有充分的理由，我想这就是其中之一。

【参考方案1】：

您可以使用submatrix views 写入连续或非连续内存：

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat foo(arma::mat A) 
    A.submat(0, 0, 1, 1).fill(1.0);
    return A;


// [[Rcpp::export]]
arma::mat bar(arma::mat A) 
    arma::uvec rows;
    rows << 0 << 2;
    arma::uvec cols;
    cols << 0 << 2;
    A.submat(rows, cols).fill(2.0);
    return A;


/*** R
A <- matrix(0, 4, 4)
foo(A)
bar(A)
*/

输出：

> A <- matrix(0, 4, 4)

> foo(A)
     [,1] [,2] [,3] [,4]
[1,]    1    1    0    0
[2,]    1    1    0    0
[3,]    0    0    0    0
[4,]    0    0    0    0

> bar(A)
     [,1] [,2] [,3] [,4]
[1,]    2    0    2    0
[2,]    0    0    0    0
[3,]    2    0    2    0
[4,]    0    0    0    0

【讨论】：

谢谢，拉尔夫。这些子矩阵视图是我一直使用的，但也许我通过做一个太小的例子来误导你。我使用我提取的子矩阵进行各种矩阵计算，所以fill 只是为了在我的示例中做一些事情。所以我确实需要子矩阵作为它自己的对象存在，我所追求的是为我的子矩阵使用与大矩阵中的原始元素相同的内存。抱歉，如果我不清楚，我会更新我的帖子以澄清。