Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to sum by all diagonals of a matrix?

Tags:

r

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
like image 597
Dambo Avatar asked Dec 08 '22 16:12

Dambo


2 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
like image 83
akrun Avatar answered Dec 30 '22 03:12

akrun


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
like image 34
tblznbits Avatar answered Dec 30 '22 03:12

tblznbits