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
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
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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With