How to sum by all diagonals of a matrix?

user_4997164 0
Viewed 711

I have this matrix:

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

And I would like to sum on each diagonal, like:

[2,1]+[1,2]
[3,1]+[2,2]+[1,3]
[4,1]+[3,2]+[2,3]+[1,4] ... etc
r
5 Answers

we can try

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
user_37322710

Another approach simply subsets the matrix, reverses it, and sums the diagonal.

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

Here are some benchmarks for the first few solutions:

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's answer is by far the fastest, but I do not believe it's providing the correct results. See here:

> 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
user_27898630

Here's an approach

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

Could be modified to incorporate both diagonals if necessary

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
user_71289340

You can use col and row with logical condition:

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
user_10301100

EDIT: Actually on a second look, only @Frank seems to have the right solution (if your remove the first element). Other solutions stop at the middle of the matrix,or give different results in the benchmark.

I tried an approach by reordering the matrix and doing row sums on triangular(ish) matrices, not as fast with @d.b 's benchmark but still pretty good. I had to fill the matrix backward to get the upper.tri and lowertri functions to behave as desired.

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))))
user_22704750
Related