如何对矩阵的所有对角线求和？答案

【问题标题】：How to sum by all diagonals of a matrix?如何对矩阵的所有对角线求和？
【发布时间】：2017-09-25 16:05:28
【问题描述】：

我有这个矩阵：

pmat <- outer(rep(1/6, 6), c(rep(1/7, 5), 2/7), `*`)
           [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
[1,] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.04761905
[2,] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.04761905
[3,] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.04761905
[4,] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.04761905
[5,] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.04761905
[6,] 0.02380952 0.02380952 0.02380952 0.02380952 0.02380952 0.04761905

我想对每个对角线求和，例如：

[2,1]+[1,2]
[3,1]+[2,2]+[1,3]
[4,1]+[3,2]+[2,3]+[1,4] ... etc

【问题讨论】：

标签： r

【解决方案1】：

我们可以试试

res <- tapply(pmat, abs(col(pmat)- row(pmat) + ncol(pmat)), FUN = sum)
unname(res[-c(1, length(res))])
#[1] 0.04761905 0.07142857 0.09523810 0.11904762 0.16666667 0.14285714 
#[7] 0.11904762 0.09523810 0.07142857

pmat[1,2] + pmat[2,1]
#[1] 0.04761905
pmat[3,1] + pmat[2,2] + pmat[1, 3]
#[1] 0.07142857

pmat[4,1] + pmat[3, 2] + pmat[2,3] + pmat[1, 4]
#[1] 0.0952381

【讨论】：

你能检查一下你的计算吗？根据我的运行，您的方法提供了不正确的结果。请参阅以下内容：m = matrix(1:36, nrow=6); akrun(m)。这提供了以下输出：[1] 17 33 54 80 111 105 94 78 57，这是不正确的。
@brittenb 我认为我采用的对角线是相反的。如果你这样做，那会很好tapply(m, apply(abs(col(m)- row(m) + ncol(m)), 1, rev), FUN = sum) 或akrun(apply(m, 1, rev))
akrun(apply(m, 1, rev)) 即all.equal(as.vector(akrun(apply(m, 1, rev))[1:5]), m0nhawk(m))# [1] TRUE

【解决方案2】：

另一种方法只是对矩阵进行子集化、反转，然后对对角线求和。

set.seed(9025)
m = matrix(rnorm(36), nrow=6)

           [,1]       [,2]       [,3]       [,4]        [,5]       [,6]
[1,]  1.0339032  1.8949990  0.6848092  0.3779458 -1.04024392 -1.6582656
[2,]  1.5605713 -1.1165541  0.8024145 -0.6841023  0.53498823  0.2518854
[3,]  1.5165112  0.6029744 -1.3962418 -0.4186442 -0.02825636  1.7304643
[4,] -0.5837422  0.7484699 -0.6918617 -0.2382265  0.67711287 -1.0709329
[5,]  1.1523611 -1.5356641  0.8393955 -0.6276084  1.59600423  0.2778807
[6,] -0.4773037  0.2414511 -1.5736829  1.5345455  0.59178567  0.5837533

sapply(2:nrow(m), function(i) {
    x = m[1:i, 1:i]
    x = apply(x, 2, rev)
    return(sum(diag(x)))
})
[1]  3.455570  1.084766  1.199592 -1.219757 -4.246751

以下是前几个解决方案的一些基准：

akrun = function(pmat) {
    res <- tapply(pmat, abs(col(pmat)- row(pmat) + ncol(pmat)), FUN = sum)
    unname(res[-c(1, length(res))])
}

m0nhawk = function(pmat) {
    vals <- sapply(3:(nrow(pmat) + 1), function(j) sum(pmat[row(pmat)+col(pmat)==j]))
    vals
}

brittenb = function(pmat) {
    sapply(2:nrow(pmat), function(i) {
        x = pmat[1:i, 1:i]
        x = apply(x, 2, rev)
        return(sum(diag(x)))
    })
}

db1 = function(pmat) {
    sapply(2:NCOL(pmat), function(i) sum(diag(t(pmat[1:i, 1:i]))))
}

db2 = function(pmat) {
    sapply(2:NCOL(pmat), function(i) sum(diag(t(pmat[1:i, 1:i]))))
}

set.seed(9025)
m = matrix(runif(1e6), nrow=1e3)
microbenchmark::microbenchmark(akrun=akrun(m), 
                               m0nhawk=m0nhawk(m), 
                               brittenb=brittenb(m),
                               db1=db1(m),
                               db2=db2(m),
                               times=1)

Unit: milliseconds
     expr         min          lq        mean      median          uq         max neval
    akrun    54.70787    54.70787    54.70787    54.70787    54.70787    54.70787     1
  m0nhawk 14946.54554 14946.54554 14946.54554 14946.54554 14946.54554 14946.54554     1
 brittenb 33241.26196 33241.26196 33241.26196 33241.26196 33241.26196 33241.26196     1
      db1  9031.17296  9031.17296  9031.17296  9031.17296  9031.17296  9031.17296     1
      db2  9194.26180  9194.26180  9194.26180  9194.26180  9194.26180  9194.26180     1

@akrun 的回答是迄今为止最快的，但我不相信它提供了正确的结果。见这里：

> head(akrun(m))
[1] 0.8088537 1.4337571 1.7651191 2.2585257 3.3970193 4.7980865
> head(m0nhawk(m))
[1] 0.05956845 2.22206820 1.87148309 2.72107233 1.76718359 3.95530571
> head(brittenb(m))
[1] 0.05956845 2.22206820 1.87148309 2.72107233 1.76718359 3.95530571

【讨论】：

@d.b 我已经添加了你的指标，但你的结果也不符合我和 m0nhawk 的结果。它们与 akrun 的不同。
您应该在评论中包含 Frank 的解决方案：here。
@d.b 将其包含在下面的答案中。

【解决方案3】：

这是一种方法

sapply(2:NCOL(pmat), function(i) sum(pmat[cbind(i:1, 1:i)]))
#[1] 0.04761905 0.07142857 0.09523810 0.11904762 0.16666667

如有必要，可以修改以合并两条对角线

setNames(object = data.frame(do.call(rbind,
            lapply(X = 2:NCOL(pmat),
                FUN = function(i) cbind(sum(pmat[cbind(i:1, 1:i)]),
                    sum(pmat[cbind(1:i, 1:i)]))))),
         nm = c("anti_diagonal", "main_diagonal"))
#  anti_diagonal main_diagonal
#1    0.04761905    0.04761905
#2    0.07142857    0.07142857
#3    0.09523810    0.09523810
#4    0.11904762    0.11904762
#5    0.16666667    0.16666667

【讨论】：

【解决方案4】：

您可以将col 和row 与逻辑条件一起使用：

vals <- sapply(3:(nrow(pmat) + 1), function(j) sum(pmat[row(pmat)+col(pmat)==j]))
# [1] 0.04761905 0.07142857 0.09523810 0.11904762 0.16666667

【讨论】：

或vapply(split(pmat, row(pmat) + col(pmat)), sum, numeric(1))
或者，只是为了保存一些字符：tapply(pmat, row(pmat) + col(pmat), sum)

【解决方案5】：

编辑：实际上再看一遍，只有@Frank 似乎有正确的解决方案（如果您删除第一个元素）。其他解决方案停在矩阵的中间，或者在基准测试中给出不同的结果。

我尝试了一种方法，通过重新排序矩阵并对三角（ish）矩阵进行行求和，虽然没有@d.b 的基准测试那么快，但仍然相当不错。我必须向后填充矩阵以使 upper.tri 和 lowertri 函数按预期运行。

n <- ncol(pmat)
m <- matrix(pmat[c(rep(seq(length(pmat)-n+2,2,1-n),n-1)+rep((n-2):0,each=(n+1)))],ncol=n+1,byrow=T)
rev(c(rowSums(m*lower.tri(m,T)),rowSums(m*upper.tri(m))))

【讨论】：