How to find the sum of all anti-diagonals?

Question

I have a matrix M:

n = 3    
x=c(0.85, 0.1, 0.05)
M <- matrix(NA, n, n); 

for(i in 1:n){
for(j in 1:n){
M[i,j] = x[i] * x[j]
}}

#       [,1]  [,2]   [,3]
# [1,] 0.7225 0.085 0.0425
# [2,] 0.0850 0.010 0.0050
# [3,] 0.0425 0.005 0.0025

I need to find the sum of all anti-diagonals include M[1,1] and M[n, n]. My attemp is

d <-matrix(c(0, 1, 2, 1, 2, 3, 2, 3, 4), n)
tapply(M, d, sum)

     0      1      2      3      4 
0.7225 0.1700 0.0950 0.0100 0.0025 

The result is correct for me.

Question. How to define the entries of matrix d? May be as function over col(M) and row(M).

Answer 1

First note that outer can produce the matrix d without explicitly listing its elements.

matrix(c(0, 1, 2, 1, 2, 3, 2, 3, 4), 3)
#>      [,1] [,2] [,3]
#> [1,]    0    1    2
#> [2,]    1    2    3
#> [3,]    2    3    4
outer(0:2, 0:2, `+`)
#>      [,1] [,2] [,3]
#> [1,]    0    1    2
#> [2,]    1    2    3
#> [3,]    2    3    4

^{Created on 2022-03-24 by the reprex package (v2.0.1)}

And use it in a function.

sumAntiDiag <- function(M){
  nr <- nrow(M)
  nc <- ncol(M)
  d <- outer(seq.int(nr), seq.int(nc), `+`)
  tapply(M, d, sum)
}

n <- 3    
x <- c(0.85, 0.1, 0.05)
M <- matrix(NA, n, n); 

for(i in 1:n){
  for(j in 1:n){
    M[i,j] = x[i] * x[j]
  }}

sumAntiDiag(M)
#>      2      3      4      5      6 
#> 0.7225 0.1700 0.0950 0.0100 0.0025

^{Created on 2022-03-24 by the reprex package (v2.0.1)}

Answer 2

As you mention in your question, row(M) and col(M) can be used, although they start rows/columns at 1 rather than zero, so you need to subtract 2 (1 for each) giving:

tapply(M, row(M) + col(M) - 2, sum)
#     0      1      2      3      4 
#0.7225 0.1700 0.0950 0.0100 0.0025

Answer 3

You could do:

sapply(seq(3), function(x) seq(3) + x - 2)
#>      [,1] [,2] [,3]
#> [1,]    0    1    2
#> [2,]    1    2    3
#> [3,]    2    3    4

Or more generally,

anti_diagonal <- function(n) sapply(seq(n), function(x) seq(n) + x - 2)

For example:

anti_diagonal(6)
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    0    1    2    3    4    5
#> [2,]    1    2    3    4    5    6
#> [3,]    2    3    4    5    6    7
#> [4,]    3    4    5    6    7    8
#> [5,]    4    5    6    7    8    9
#> [6,]    5    6    7    8    9   10

Answer 4

Try the code below by defining a function f using embed from base R, i.e.,

f <- function(n) embed(seq(2 * n - 1) - 1, n)[, n:1]

such that

> f(3)
     [,1] [,2] [,3]
[1,]    0    1    2
[2,]    1    2    3
[3,]    2    3    4

> f(4)
     [,1] [,2] [,3] [,4]
[1,]    0    1    2    3
[2,]    1    2    3    4
[3,]    2    3    4    5
[4,]    3    4    5    6

> f(5)
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    1    2    3    4
[2,]    1    2    3    4    5
[3,]    2    3    4    5    6
[4,]    3    4    5    6    7
[5,]    4    5    6    7    8

Answer 5

You can use sequence:

function(n) matrix(sequence(rep(n, n), seq(n) - 1), nrow = n)

output

f <- function(n) matrix(sequence(rep(n, n), seq(n) - 1), nrow = n)
f(3)
     [,1] [,2] [,3]
[1,]    0    1    2
[2,]    1    2    3
[3,]    2    3    4

f(5)
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    1    2    3    4
[2,]    1    2    3    4    5
[3,]    2    3    4    5    6
[4,]    3    4    5    6    7
[5,]    4    5    6    7    8

Answer 6

Using indexing instead of tapply will speed things up a bit. Or Rcpp:

sumdiags <- function(mat, minor = TRUE) {
  m <- ncol(mat)
  
  if (minor) {
    n <- nrow(mat)
    lens <- c(pmin(1:n, m), pmin((m - 1L):1, n))
    c(mat[1], diff(cumsum(mat[sequence(lens, c(1:n, seq(2L*n, by = n, length.out = m - 1L)), n - 1L)])[cumsum(lens)]))
  } else {
    Recall(mat[,m:1])
  }
}

# compare to tapply solution
sumdiags2 <- function(mat, minor = TRUE) {
  if (minor) {
    as.numeric(tapply(mat, row(mat) + col(mat), sum))
  } else {
    Recall(mat[,ncol(mat):1])
  }
}

# or Rcpp
Rcpp::cppFunction('NumericVector sumdiagsRcpp(const NumericMatrix& mat) {
  const int n = mat.nrow();
  const int m = mat.ncol();
  NumericVector x (n + m - 1);

  for(int row = 0; row < n; row++) {
    for(int col = 0; col < m; col++) {
      x[row + col] += mat(row, col);
    }
  }
  return x;
}')

# OP data
x <- c(0.85, 0.1, 0.05)
m <- outer(x, x)
sumdiags(m)
#> [1] 0.7225 0.1700 0.0950 0.0100 0.0025
sumdiags2(m)
#> [1] 0.7225 0.1700 0.0950 0.0100 0.0025
sumdiagsRcpp(m)
#> [1] 0.7225 0.1700 0.0950 0.0100 0.0025

# bigger matrix for benchmarking
m <- matrix(runif(1e6), 1e3)
microbenchmark::microbenchmark(sumdiags = sumdiags(m),
                               sumdiags2 = sumdiags2(m),
                               sumdiagsRcpp = sumdiagsRcpp(m),
                               check = "equal")
#> Unit: milliseconds
#>         expr       min        lq      mean    median      uq        max neval
#>     sumdiags  9.985302 10.266350 13.686723 10.803401 17.5274  22.387601   100
#>    sumdiags2 55.790402 65.140051 78.763478 67.120051 70.4165 183.936801   100
#> sumdiagsRcpp  2.192201  2.378651  2.599326  2.631751  2.7050   4.038301   100