将矩阵下对角线转换为列 r 的有效算法

Question

我有一个非方阵，需要对其下对角线进行一些计算。我发现最好的方法是将下对角线转换为列/行并使用 cumprod 之类的函数。现在我使用定义如下的 for 循环和 exdiag：

exdiag <- function(mat, off=0) {mat[row(mat) == col(mat)+off]}

但是它并不是真正有效的。你知道任何其他算法可以实现这种结果吗？

一个小例子来说明我在做什么：

exdiag <- function(mat, off=0) {mat[row(mat) == col(mat)+off]}
mat <- matrix(1:72, nrow = 12, ncol = 6)
newmat <- matrix(nrow=11, ncol=6)
for (i in 1:11){
   newmat[i,] <- c(cumprod(exdiag(mat,i)),rep(0,max(6-12+i,0)))
}

最好的问候， 阿图尔

Answer 1

从非方阵中获取所有可能对角线的最快但迄今为止最神秘的解决方案是将矩阵视为向量并简单地构造一个id 向量以供选择。最后，您可以根据需要将其转换回矩阵。

以下函数可以做到这一点：

exdiag <- function(mat){

  NR <- nrow(mat)
  NC <- ncol(mat)
  smalldim <- min(NC,NR)
  if(NC > NR){
    id <- seq_len(NR) + 
      seq.int(0,NR-1)*NR +
      rep(seq.int(1,NC - 1), each = NR)*NR


  } else if(NC < NR){
    id <- seq_len(NC) + 
      seq.int(0,NC-1)*NR +
      rep(seq.int(1,NR - 1), each = NC)

  } else {
    return(diag(mat))
  }
  out <- matrix(mat[id],nrow = smalldim)
  id <- (ncol(out) + 1 - row(out)) - col(out) < 0
  out[id] <- NA
  return(out)
}

请记住，您必须考虑矩阵的形成方式。

在这两种情况下，我都遵循相同的逻辑：

• 首先构造一个序列，指示沿最小维度的位置
• 在此序列中，添加 0、1、2、... 乘以行长。

这将创建第一个对角线。完成此操作后，您只需添加一个序列，将整个前一个序列移动 1（向下或向右），直到到达矩阵的末尾。要右移，我需要将此序列乘以行数。

最后，您可以使用这些索引从 mat 中选择正确的位置，并将所有这些作为矩阵返回。由于此代码的矢量化特性，您必须检查最后的下对角线是否正确。这些包含的元素比第一个少，因此您必须将不属于该次对角线的值替换为 NA。同样在这里，您可以简单地使用索引技巧。

您可以按如下方式使用它：

> diag1 <- exdiag(amatrix)
> diag2 <- exdiag(t(amatrix))
> identical(diag1, diag2)
[1] TRUE

为了得出你的结果

amatrix <- matrix(1:72, ncol = 6)
diag1 <- exdiag(amatrix)
res <- apply(diag1,2,cumprod)
res[is.na(res)] <- 0
t(res)

Answer 2

您可以修改 diag() 函数。

exdiag <- function(mat, off=0) {mat[row(mat) == col(mat)+off]}

exdiag2 <- function(matrix, off){diag(matrix[-1:-off,])}

速度测试：

mat = diag(10, 10000,10000)

off = 4

> system.time(exdiag(mat,4))

  user  system elapsed 
  7.083   2.973  10.054 

> system.time(exdiag2(mat,4))

   user  system elapsed 
  5.370   0.155   5.524 

> system.time(diag(mat))

   user  system elapsed 
  0.002   0.000   0.002

看起来矩阵的子集需要很多时间，但它仍然比您的实现更好。可能还有很多其他子集方法，它们优于我的解决方案。 :)